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ABSTRACT 

Wc perform a suite of two- and three-dimensional magnctohydrodynamic (MHD) sim- 
ulations with the Athena code of the Kelvin-Helmholtz instability (KHI) in the sub- 
sonic, weak magnetic field limit and study the non-linear development leading to 
magnetic field growth, saturation, and subsequent turbulent decay. Two-dimensional 
KHI simulations do not converge beyond the linear growth stage, while their three- 
dimensional counterparts converge across all stages of evolution. Appealing to spectral 
energy transfer function analysis, wc quantify energy transfer on a scale-by-scale basis 
and identify the physical mechanisms responsible for energy exchange. At late times 
when the fluid is in a state of MHD turbulence, magnetic tension mediates the dom- 
inant mode of energy injection into the magnetic reservoir, whereby turbulent fluid 
motions twist and stretch the magnetic field lines. This generated magnetic energy 
turbulently cascades to smaller scales, while being exchanged backwards and forwards 
with the kinetic energy reservoir, until finally being dissipated. Extending the ideal 
MHD treatment to include explicit viscous and resistive terms, we demonstrate that 
physical dissipation terms, as opposed to numerical effects, determine the dissipation 
scale. For ideal MHD simulations, the dissipation scale shifts to progressively smaller 
scales as numerical resolution is increased; however, the inclusion of physical dissipa- 
tion reduces the importance of numerical dissipation and acts to move the dissipation 
scale to larger spatial scales. This is because physical dissipation acts across all scales, 
while numerical dissipation operates preferentially on the (generally small) dissipa- 
tion scale imposed by the grid resolution; therefore, incorporating explicit dissipation 
pushes the dissipation scale to larger scales than if the dissipation were entirely nu- 
merical. For scales larger than the dissipation scale, we show that the physics of energy 
transfer in decaying MHD turbulence are robust to numerical effects. 
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1 INTRODUCTION 

Shearing flows potentially susceptible to the Kelvin- 
Helmholtz instability (KHI) are ubiquitous in astrophysi- 
cal systems ranging from black hole jets to protoplanetary 
accretion discs. Given the tendency of the KHI to disrupt 
shear layers and generate turbulence, it is essential that we 
understand in which astrophysical contexts the instability is 
expected to occur, how quickly it grows, and how it even- 
tually saturates. These issues gain an added significance in 
the presence of magnetic fields, where the KHI can work 
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to enhance weak fields through a combination of twisting, 
shearing, and compression until the magnetic fields become 
strong enough to take a more active role in their evolution. 
If a magnetic field component oriented along the direction 
of the flow is sufficiently strong, magnetic tension resists 
mixing across the shear layer, which in turn limits further 
development of the KHI. 

Whether an astrophysical flow is unstable to the KHI 
will depend upon multiple parameters. Highly relativis- 
tic flows, such as those inferred to exist in gamma-ray 
bursts and observed in some active galactic nuclei and 
blazars, ar e not l i kely to be strongly affected by the 
KHI (e.g., lHardeej l2007l h Systems in which portions of 
the flow are significantly subluminal, however, are more 
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likely to become unstable, and the KHI has been invoked 
to explain the knots in M87 (|Bicknell fc Begelmanl Il996l ; 
lHardee fc Eilekl l201ll) and similar features in other jet- 
ted s ystems IjHardee. Walker fc Gom ez 20051; IPerucho et al.l 
120061 ). Shearing subsonic and sub-Alfvenic flows will gener- 
ally be unstable to the KHI, with growth rates and satu- 
ration behaviours that depend upon the exact physical cir- 
cumstances under consideration. 

Much of the seminal work on understanding the linear 
growth of the KHI was undertaken in generalised studies 
of terrestrial fluid mechan i cs, such as t hos e considered by 
Lamb dl932h. iMilesl dl957l . ll95Sl Il959al lbh. IChandrasekharl 
(|l96lh . and lFeier fc Miles! \ 19631 ). The first properly astro- 
physical application of KHI ph ysics appear s to have been to 
the E arth's magnetosphere by ISenl ( 19651 ) and ISouthwoodl 
() 19681 ') . with a general review of various combina tions of 
compr e ssibilit y and magnetisation provided by iGerwinl 
l|l968l ). lBlak3 (|l972h considered KHI growth in radio lobes, 
followed by treatments of KH I in relativistic beams b y 
iBlandford fc Pringld | 19761 ) a ndlTurland fc Scheuerl (| 19761). 



Since t hat time, many grou p s ([Ferrari. Trussoni fc Zaninettil 
19781; iMiura fc Pritchettl Il982l; iBodo fc Ferraril Il982l 



Birkinshawlll984l: iFiedler fc Jonedll984l; lHardee fc Normanl 



19881; iBirkinshawl Il99ll; IAppI fc Camenziiidl 1 19921 ; 
Hardee fc Stonel Il997l ; lHarded |2000| . 120071 ) ""have pro- 



vided refined analytic treatments of KHI development 
in jets, including considerations such as geometry, shear 
profile, magnetisation, and radiation. 

While these studies helped to pin down the linear evolu- 
tion of the KHI, numerical methods have proven essential for 
understanding the non-linear development of the instability. 
In two spatial dimensions (2D), hydrodynamics si mula- 
tions of sh e aring flows by iNorman fc Harcleel (1 19881 ) and 



IBodo et all (jl994ll995l ) provided some of the first glimpses 
into the non-linear evolution of the KHI, followed by exten- 
sio ns into 2D/ax i symm et ric magnetoh y drodyn a mics (MHD) 
by iFrank et all (1 19961) I Jones et~aT1 (Il997al). Ijeong et al.l 
(|2000l ). IStone fc Harded (|2000l ). and IPalotti et all (|2008l ). 
Full three-dimensional (3D) nume rical explorations 



of th e KHI were conducted b y Norman fc Balsaral 



(|l993l ) and iBassett fc Woodward! Jl995l) in th e h; 



drodynamic limit and 



Hardee. Clarke 



m tn e ny - 
Rosenl dl997l ). 



iRvu. Jones fc Frank! d2000h. and lHardee fc Rosen! I2OO2I) 



full MHD. Mo re recen tl y iMartf. Perucho fc Hanasj 



Perucho et al 



Perucho et al 



2004) 



(120041 ); IPerucho. Martf fc Hanasj (|2004| ). 
20061 ). and iRadice fc Rezzollal (|2012l ) have 
explored the KHI in rela t ivistic hyd rodynamics, while 
IZhang, MacFadven fc Wand (|2009l ) and iBeckwith fc Stond 



( 201 ll ) have discussed KHI development in full relativistic 
MHD. 

The majority of these numerical studies analysed 
KHI development through the measurement of instabil- 
ity growth rates, saturation behaviours, and/or morpho- 
logical consequences of instability. In this study, we will 
provide a novel look at the development of the KHI 
by employing sp e ctral energy transfer function analysis 
dKraichnanl 1 19671 : IPietarila Graham. Cameron fc Schiisslerl 
120101 ) to determine precisely how energy is transferred in 
the KHI to specific spatial scales as a function of en- 
ergy type (kinetic, magnetic, internal) and transfer force 
(compression, advection, magnetic tension, magnetic pres- 
sure). This approach provides us with insight into the de- 



tails of the KHI physics that are not otherwise accessible 
and allows us to determine how integrated flow proper- 
ties and morphology reflect the scale-dependent processes 
we identify. Transfer function analysis has far-reaching 
astrophysical app lications including studying the turbu- 
lent s olar dynamo ijPietarila Graham. Cameron fc Schusslerl 
120101 ). understanding the detailed physics of accretion 
disc tu rbulence arising from the m a gnetorotational insta- 
bility jFromang fc Papaloizoul 120071 : iFromang et al.l 120071 : 
bimon, Hawlev fc Beckwithll2009l ). and quantifying the dy- 



namical importance of "mesoscale" magnetic structures that 
arise in l ocal studies of accretion discs with large spatial 
domains dSimon. Beckwith fc Armitagd I2012T ) . This treat- 
ment also provides a scale-by-scale look into the properties 
of energy dissipation distributions, which are of particular 
relevance in considerations of the sensitivity of emergent 
spectra from black hole systems to the accretion disc atmo- 
sphere and coronal structure (e.g.. iMerloni. Fabian fc Ross 
20001; lTurnerll2004l ; iHirose. Krolik fc Stonel2006l ; lBlaes et al 



20061 ) . Additionally, we will discuss how computational is- 



sues, such as numerical convergence, the effects of domain 
size, and numerical dissipation behaviours, can be under- 
stood and evaluated in terms of KHI development. 

While the KHI has important applications to subsonic, 
transonic, supersonic, and relativistic astrophysics, we fo- 
cus here on understanding the non-linear development and 
spectral structure of the KHI in the subsonic, weakly mag- 
netised limit. The motivations for this choice are both sim- 
plicity and applicability. We wish to apply comprehensive 
analysis tools to study the development of the KHI for a sim- 
plified case without the complications of additional physics 
like shock formation or the exchange between different fluid 
instabilities such as the family of current-driven instabilities 
(|Begelmanlll998l ). Particular attention is given to properly 
constructing an initial setup for the simulations and provid- 
ing convincing evidence that the late-stage development is 
physical, rather than numerical, in origin. We aim for our 
numerical study of the KHI to be relevant and extendable 
to a broad range of astrophysical applications, such as the 
interplay between the KHI and CDI in jets, the nature of 
MHD turbulence arising from the KHI, and dissipation pro- 
files in accretion discs. 

This paper is organised as follows. Descriptions of the 
Athena MHD code and the KHI problem setup are provided 
in ^2] The methodology behind the spectral energy transfer 
function analysis we adopt is given in fj3j The convergence 
of 3D KHI simulations is discussed in §Z\ along with the 
inadequacy of 2D simulations. We next describe the evolu- 
tion of the KHI simulations with a focus on the late-stage 
turbulent decay and the importance of energy transfer in- 
volving the magnetic energy reservoir in ij5] The inclusion of 
physical dissipation is explored in [J6] Finally, fJ7| is devoted 
to a summary and discussion of this work, followed by our 
conclusions in 351 



2 NUMERICAL DETAILS 



We solve the equ ations of MHD using the Athena code 
(|Stone et al.ll200a ). a second-order accurate Godunov flux- 
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conservative codtQ. Athena is an Eulerian code that solves 
the equations of compressible, adiabatic, inviscid, ideal 
MHD in conservative form, 



2lL 

dt 

d(pv) 
dt 
dE 
~dt 

m 



= - V • (pv) 



V- 



pvv BB 



(E + P+^B 2 ^j v B(B v] 



— = V x (v x B) . 



(1) 

(2) 
(3) 
(4) 



The notation is of familiar form, where p is the density, v is 
the fluid velocity, P is the gas pressure, B is the magnetic 
field, and E is the total energy density defined by, 



E = e+-pv +2 B > 



(5) 



where e = P/ (7 — 1) is the internal energy density for an 
ideal gas and 7 is the adiabatic index. I is the identity matrix 
operating on the total pressure, P+B 2 /2. In the adopted no- 
tation, the magnetic field absorbs a factor of ^Jp/in, where 
fi = 1 is the assumed magnetic permeability. The MHD 
equations Q~HJ] are conservation equations describing, in or- 
der, the conservation of mass, momentum, total energy, and 
magnetic flux. Equations [T}[3] have the generic form of any 
conservation equation, where the time derivative of a con- 
served quantity is equated to the divergence of a flux, in the 
absence of any source/sink terms. Dissipation terms such 
as viscosity, resistivity, and conduction are neglected in our 
treatment of ideal MHD. We investigate the addition of ex- 
plicit dissipation terms and their affect on the KHI evolution 
in SHU 

The b asic algorithms implemented i n Athena are de- 
scribed by iGardiner fc Stone! (|2005l . 12008ft with further de- 
tai ls (impleme n tation and multi-dimensional tests) given 
in IStone et al.l (|200ct ). Specifically, we use the dimension- 
ally unsplit Corner Tran s port Upwind (CTU) integrator 
described by IStone et alj (|2008l ) combined with the con - 
strained transport (CT) method of lEvans fc Hawlevl ll 19881 ) 
to maintain the divergence-less nature of the magnetic field 
in multi-dimensions. Athena implements a variety of op- 
tions for spatial reconstruction and solution of the Riemann 
problem. In this work, we use third-order spatial recon- 
struction in characteristic variables and the HLLC/HLLD 
Riemann solvers for hydrodynamic/MHD simulations. We 
avoid choosing the HLLE solver due to its highly diffusive 
behaviour (for further information, see Appendix[S]|. In this 
work, we make extensive use of the conservation properties 
of Athena to examine exchange of energy between kinetic, 
magnetic, and internal energy reservoirs. 



2.1 Problem Setup 

The initial problem setup for numerical simulations of the 
KHI with Athena is shown schematically in Figure [T] We 
consider a 3D Cartesian grid centered on (x, y, z) = (0, 0, 0) 



with dimensions L x 



1 and periodic 



1 The Athena code and a repository of test problems are main- 
tained online at https://trac.princcton.edu/Athena/. 
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Figure 1. Schematic of the KHI problem setup for the Athena 
simulations. Each side of the computational box has length L, 
with the origin at the center, (x, y, z) = (0, 0, 0). Periodic bound- 
ary conditions are adopted in all directions. Counter-streaming 
flows are initiated in the {/-direction, each with speed Uo in the 
laboratory frame. The box is filled with uniform pressure, Po, and 
a uniform magnetic field of strength Bo in the {/-direction. The 
fluid densities are pi and P2 in Region 1 and Region 2, respec- 
tively. Region 1 is bounded in the 2-direction by — zo < z < zo and 
Region 2 corresponds to z > |zo|i where the shear interfaces are 
located at rbzo- Although not represented in the schematic, the 
density and velocity profiles across each shear layer are matched 
continuously by a hyperbolic tangent function; thus, permitting 
the interfaces to be well-resolved. 



boundary conditions enforced in all directions. Counter- 
streaming flows are set up along the y-direction according 
to the velocity profile, 

J i/otaAffa), \z\>zo 
[ -[/ tanh(^aJ^J , \z\ < zo 

where 2Uo = 1 is the magnitude of the relative shear veloc- 
ity, zo = 0.25L specifies the location of the shear interfaces, 
and a — 0.01L is a parameter describing the thickness of 
the shear layer. A continuous velocity profile is constructed 
across the shear layers, rather than a discontinuous inter- 
face, to ensure that truncation error resulting from numer- 
ical diffusion of unresolved modes (i.e., short wavelength, 
large wavenumber) does not dominate the solutions (see 
Appendix . The linear growth rate of the KHI is pro- 
portional to the wavenumber; thus, an under-resolved shear 
layer will evolve unphysically into the non-linear regime. The 
hyperbolic tangent profile we adopt provides a sharp, yet 
smooth, transition while also introducing the length scale, 
a, to an otherwise scale-free problem. For a given grid res- 
olution, N = N x = N y = N z , the shear layer is resolved 
by AaN/L grid zones, which amounts to ~ 20 resolution 
elements across the interface for the = 512 3D MHD 
simulation, which is the fiducial run for the majority of the 
analysis. The wavenumber corresponding to the full width of 
the shear layer is fc s h = 27r/(4a) ~ 157; therefore, to resolve 
modes that grow on the same spatial scale of the shear layer 
or smaller, the simulation resolution must be adequate, such 
that the Nyquist wavenumber, kn y = (N/2)(2n / L), exceeds 
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Region 


Vy 


P 


P 


c 3 


J( 


Bo 


A) 


OA 




1 


-0.5 


1 


2 


0.91 


0.39 


0.02 


5000 


0.014 


35.4 


2 


0.5 


1 


1 


1.29 


0.55 


0.02 


5000 


0.020 


25.0 



Table 1. Initial conditions of the KHI problem setup for Athena 
simulations. From left to right the columns are the Region number 
(see Figure[T]l, shear flow velocity, gas pressure, gas density, sound 
speed, Mach number, strength of the magnetic field aligned with 
the shear flow, gas-to-magnetic pressures ratio, Alfvcn speed, and 
Alfven Mach number. Where units apply but are left unspecified, 
these are arbitrary code units. 



The initial density profile is described by. 
p(z) = 



2 Vp 2 



20 



- 1 



where p\ — 2 is the density of the central fluid slab and 
p2 = 1 is the density of the surrounding fluid. The con- 
tact discontinuity is smeared by the same hyperbolic tangent 
function applied to the velocity profile to ensure a resolved 
solution. Initially, the fluid slabs are in gas pressure equilib- 
rium with Pq = 1 and adiabatic index, 7 = 5/3. A uniform 
magnetic field, Bo = Boy, is aligned parallel to the shear 
flow with initial strength, Bo = 0.02, corresponding to the 
weakly non-linear regime, meaning that the magnetic field is 
weak and the flow is not linearly stable. In this regime, the 
instability is essentially hydrodynamic early on, then enters 
the non-linear regime where secondary instabilities break up 
large-scale structures and magnetic energy is amplified due 
to twisting/stretching of magnetic field. After saturation, 
the flow enters a state of decaying MHD turbulence (for a 
discussio n of different stability regimes of the magnetise d 
KHI, see iRvu. Jones fc Frank|[200oTlBatv fc Keppenj|2002h . 

In order to provoke the onset of the KHI, at time t = 
we impose a small-amplitude, single mode perturbation to 
the z-component of velocity of the form, 



v' z (x, y, z) = «°sin (k x x) sin (k y y) e 



-(z + z ) 2 /a 



(8) 



where v® — O.OlcTo, k x — k y = 2ir/L, and a = 0.1L de- 
scribes the decaying behaviour of the perturbation along the 
z-direction. A full perturbation wavelength fits within the 
x and y computational box boundaries. Modes with wave- 
lengths larger than the box scale, L, are not captured. 

Table [T] summarises the set of initial parameters corre- 
sponding to each region defined in Figure [1] All simulations 
used the foregoing setup and parameter choices, unless spec- 
ified otherwise. A parameter survey is beyond our scope and 
does not address our motivating intention to study in detail 
the development and energetics of the KHI starting from a 
properly constructed initial configuration. Table [2] lists the 
suite of simulations presented in this work. 



ID a 


N 


^start 




Bo 




V/V&d 


3M1024 


1024 


to 
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0.02 








3M512 


512 


to 


1 


0.02 








3M256 


256 


to 


1 


0.02 








3M128 


128 


to 


1 


0.02 








2M16384 


16,384 


to 


1 


0.02 








2M8192 


8192 


to 


1 


0.02 








2M4096 


4096 


to 


1 


0.02 








2M2048 


2048 


to 


1 


0.02 








2M1024 


1024 


to 


1 


0.02 








2M512 


512 


to 


1 


0.02 








2M256 


256 


to 


1 


0.02 








2M128 


128 


to 


1 


0.02 








3H512 


512 


to 


1 











2H8192 


8192 


to 


1 











3M512^ 


512 


to 


1 


0.02 


2 


2 


3M512^ 


512 


to 


1 


0.02 


2 


1 


3M512^ 


512 


to 


1 


0.02 


1 


2 


3M512}^ 


512 


to 


1 


0.02 


1 


1 


3M512D^ 


512 


^pcak 


1 


0.02 


2 


2 


3M512D}^ 


512 


^pcak 


1 


0.02 


1 


1 


3M256D}^ 


256 


^pcak 


1 


0.02 


1 


1 


3M128D}^ 


128 


^pcak 


1 


0.02 


1 


1 


3M512D"£j 

T//2 


512 


^pcak 


1 


0.02 


1/2 


1/2 


3M512D"^ 

77/4 


512 


^pcak 


1 


0.02 


1/4 


1/4 


3M512D"^ 


512 


^pcak 


1 


0.02 


1/8 


1/8 


3M512J 6 


512 


to 


1 


0.02 








3M512z2 


512 c 


to 


2 


0.02 








3M512z4 


512 c 


to 


1 


0.02 









Table 2. Table of KHI simulations referred to in our study. From 
left to right the columns are the simulation identification tag, grid 
resolution in each dimension, time at which the simulation was 
started from, size of the z-domain (in code units), initial magnetic 
field strength (in code units), kinematic viscosity coefficient rela- 
tive to the fiducial value, and Ohmic resistivity coefficient relative 
to the fiducial value. "The ID tag generally follows a straightfor- 
ward, three-part naming convention. The first number indicates 
the dimensionality (i.e., 2D or 3D), the letter denotes the gas 
dynamics used (i.e., M for MHD or H for hydrodynamics), and 
the trailing number specifics the grid resolution in each direction. 
''The shear layer in this simulation was discontinuous, correspond- 
ing to the width parameter, a = 0. All other simulations adopt 
a = 0.01 (in code units). Simulations 3M512z2 and 3M512z4 
have N z = 1024 and N z = 2048, respectively. 



3 SPECTRAL ANALYSIS 

Throughout this work, we use energy power spectra to ex- 
amine at which scales the majority of the magnetic energy is 
generated and how the spectral shape of the different energy 
reservoirs (i.e., kinetic, magnetic, and internal) evolve. The 
kinetic, magnetic, and internal energy power spectra, also 



referred to as spectral energy densities, are defined as, 



E K (k) = ~ (v(fc) • [pvf (fc)) 
E M (k) = i (B(fc).B*(fe) 



Ei(k) = 



P(k) 



r 



(0) 
(10) 

(ii) 
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Figure 2. Schematic representation of how differential spectral 
energy density is spherically integrated over concentric fc-shells of 
constant thickness, Afc. For a given spectral energy density, E(k), 
the differential spectral energy density, AE(k i+ i/2), is computed 
as the sum-total energy contained within a shell of inner bound- 
ary ki and outer boundary fci+i, where Afc = k%-\-i — ki is con- 
stant. To avoid double-counting, the spectral energy contained 
at the outer boundary of a shell (i.e., £(fcj_|_i) for AE^^)) 
is omitted, while the energy at the inner boundary is taken to 
be inclusive. Spectral energy is integrated over shells of thick- 
ness AfcL/(27r) = 1 from fc m j n L/(27r) = to fc max L/(27r) = 
fcN y L/(27r) = N/2. Plots of energy power spectra are showing 
the differential spectral energy density contained within a fc-shell, 
dE(k)/dk = AE(k)/Ak. 



where k = |k| = yjk% + ky + k\ , an asterisk superscript 
denotes a complex conjugate, and F(k) indicates the Fourier 
transform of the quantity /(x), 



F(k) 



/(x)e- !k ' x d 3 x. 



(12) 



No normalisation is performed, as the magnitude of energy 
transfer is of interest, rather than merely the spectral shape. 
To improve statistics and aid in interpretation, the energy 
power spectra are integrated over concentric spherical shells 
of thickness, AfcL/(27r) = 1, as shown schematically in Fig- 
ure This yields the differential power contained within a 
shell, 



dE(k) _ AE(k) 



(13) 



dk Ak ' 

centered on half-integer values of wavenumber k. 

Sometimes we will have cause to perform a shell 
averaged normalisation of spectral energy density (see 
according to, 



dE((k)) 
dk 



fc^dfc 
dk 



dE{k) 
dk ' 



(14) 



where the angled bracket convention, dE ({k)) /dk, indicates 
that a shell average was performed on the spectral energy 
density. 

Power spectra describe the distribution of energy across 
spatial scales; however, such distributions provide no clear 
way of determining how the energy transfers between 
scales and different forms (e.g., between magnetic and ki- 
netic energies). Transfer function analysis, first introduced 



by iKraichnanl Jl967h . allows for the scale-by-scale quan- 
tification of energy transfer between reservoirs and iden- 
tification of the mechanism responsible for the energy 
exchange. The mechanics of deriving the transfer func- 
tions are given in Appendix [B] with an outline of the 
derivations and an explanation of notation given here. In 
this, we closely follow the approach and inter operation of 
IPietarila Graham. Cameron fc Schiisslerl |2010l ). who rigor- 
ously developed transfer analysis for compressible MHD in 
the context of the small-scale solar turbulent dynamo. We 
specialise to the case of the KHI by decomposing the velocity 
field into cont ributions from the shearing flow and turbulent 
fluctuations ([Fromang fc Papaloizoul [20071 ; iFromang et al.l 
120071 ; ISimon, Hawlev fc Beckwithll2009f ). 



V sh + V t , 



(15) 



where v t is the turbulent velocity field and v s h is the back- 
ground flow field, defined as, 



v sh = v sh (z) y = 



L X Ly 



v v (x,y,z)dxdy. 



(16) 



Inserting the decomposed velocity field into the momen- 
tum, energy, and induction equations, taking the Fourier 
transform, and performing the appropriate dot product (see 
Appendix [B]) . the complete transfer function equations for 
kinetic, magnetic, and internal energies are, 



dE K (k) 
dt 



dE M (k) 
dt 



dEi(k) 
dt 



= TlKc(fc) + SlKc(fc) + T K KA(fc) + X K KA(k) + 

TbktW + SbktM + T B Kp{k) + S B Kp(fc) + 
TKKc(fc) + SKKc(fc) + XnKc(k) + D K (k) (17) 

= TBBA(fc) + 5*BBA (fc) + TkBT(fc) + S'KBT(fc) + 



IkBp(fe) + D M (k) 



(18) 



TkiA(fc) + SkiaM + Trig (fe) + Di(k), (19) 



where the notation is described below. These are time evo- 
lution equations of spectral energy densities. Fourier trans- 
forms are computed according to Equation [12] using a fast 
Fourier transform algorithm. The shear layer is not driven 
and is continually decaying; thus, the energy densities in the 
saturated state are not in a steady state and time derivatives 

are calculated exp licitly. 

F ollowing IPietarila Graham. Cameron fc Schiisslerl 
(|2010t ). we interpret the transfer function T, S, JfxYp(fc) 
as measuring the net energy transfer rate from all scales 
of reservoir X to scale fc of reservoir Y, where the energy 
exchange is via the force F. The net energy transfer from 
reservoir X into reservoir Y at scale fc is positive (nega- 
tive) for T, S, XxYp(fc) > (< 0). The available energy 
reservoirs are kinetic (K), magnetic (M), and internal 
(I). The mediating forces (F) depend on the exact form 
of each transfer function, but in general these forces are 
compressive motions (C), advection (A), magnetic tension 
(T), and magnetic pressure (P). Energy transfer due to 
purely turbulent motions, v t , is denoted by Txyf(&); the 
background shear flow, v s h, by SxYF(fc); or some hybrid 
cross term involving both v t and v s h by Xxyp(A;). Finally, 
the terms D^{k), Dyi{k), and Di(fc) in Equations [TTHT^l are 
simply the residuals of the time derivative of spectral energy 
density and the sum of all transfer functions, resulting in 
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Figure 3. Volume-averaged rms velocity component transverse 
to the shear layers for the resolved shearing runs 3M1024 (dot- 
ted line), 3M512 (solid line), 3M256 (dashed line), and 3M128 
(dash- dot line). Convergence is demonstrated in the linear growth 
regime for all resolutions considered here, but (n 2 ) 1 / 2 cannot be 
used as a diagnostic of convergence in the decaying regime. Here, 
and in all subsequent figures, time is paramctcriscd in units of 
the linear growth e-folding time, t, as computed from the fiducial 
3M512 simulation. 



a mea sure of numerical dissipation rate as a function of 
scale dFromang fc Papaloizou _ 20071 ; iFromang et al.l 120071 ; 
ISimon. Hawlev fe Beckwithll2009l ) 

All transfer functions are spherically integrated over 
shells of constant thickness AkL/(2ir) = 1 and then plotted 
as k- (dTxYF(k) / dk) vs. log(fc) so that a peak in the spectrum 
corresponds to the wavenumb er containing the most power 
(|Zdziarski fc Gierliriskil |2004| ). We choose to time-average 
the transfer functions over the same intervals shown in the 
energy power spectral analysis of Figure [TO] in order to im- 
prove statistics across all k and make robust statements re- 
garding energy exchange during different stages of the KHI 
evolution. 



4 CONVERGENCE 

Convergence of the linear growth stage of the KHI can 
be assessed through the volume-averaged root mean square 
(rms) velocity tr ansverse to the shear layer, (v 2 ) 1 ^ 2 (e.g., 
iFrank et alj|l996h . In this work, angled brackets surround- 
ing a quantity denote volume averages, where the volume- 
average of quantity Q(x, y, z) is given by, 



(Q) 



L X LyL Z 



1 (x, y, z) dxdydz. 



(20) 



Figure |3] shows {v 2 ) 1/2 for the runs 3M128, 3M256, 3M512, 
and 3M1024, where the initial value is dictated by the 
perturbation of the equilibrium configuration. The linear 
growth phase of the KHI is identified as the exponentially in- 
creasing portion of (wf} 1//2 . We choose to parameterise time 
in terms of the linear growth e-folding time, r ~ 0.16i, eval- 
uated over the exponentially increasing portion of (v 2 ) 1 ^ 2 
according to, (v 2 ) 1 ^ 2 = Ae i//r , where A is the initial rms 
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Figure 4. Time evolution of 3D MHD simulations showing the 
amplification, saturation, and decay of the volume-averaged mag- 
netic field magnitude, (B), relative to the initial volume-averaged 
magnetic field magnitude, (-Bo)- The 3D simulations shown in 
this convergence study are 3M128 (dash-dot line), 3M256 (dashed 
line), 3M512 (solid line), and 3Mi024 (dotted line). The times 
marked by vertical lines correspond to the 3M512 simulation and 
are defined in the text (see jj5}. Convergence is demonstrated at 
a resolution, N = 512, for the linear growth and non-linear decay 
regimes, with only modest differences in the saturation amplitude. 
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Figure 5. Time evolution of 2D MHD simulations showing the 
amplification, saturation, and decay of the volume-averaged mag- 
netic field magnitude, (B), relative to the initial volume-averaged 
magnetic field magnitude, (-Bo). The 2D simulations shown in this 
convergence study are 2M128 (red line), 2M256 (orange line), 
2M512 (green line), 2M1024 (cyan line), 2M2048 (blue line), 
2M4096 (violet line), 2M8192 (black line), and 2M16384 (pur- 
ple line). The times marked by vertical lines correspond to the 
2M8192 simulation and are defined in the text (see jjjj. While the 
linear regime is well-converged at a resolution, N = 512, neither 
the non-linear saturated state nor the late-time decay show in- 
dications of convergence, even for the very high resolution case, 
N = 16,384. 
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Figure 6. Time evolutions of the volume-averaged quantities: 
internal energy, (Ei) (top panel), kinetic energy, (Bk) (middle 
panel), and magnetic energy, (Em) (bottom panel). All volume- 
averaged energies are shown relative to the total volume-averaged 
energy in the computational box, (-Etot). Shown are results for 
KHI simulations 3M1024 (black dotted lines), 3M512 (black solid 
lines), 3H512 (red solid lines), 2M8192 (black dashed lines), and 
2H8192 (red dashed lines). Magnetic energy is more efficiently 
generated from the available kinetic energy and less efficiently 
dissipated for the 2D KHI. 



transverse velocity and t is time in code units. Although 
r is decreasingly relevant as the flow becomes turbulent, 
it is physically motivated and well-defined during the lin- 
ear growth. Figure [3] shows that the linear growth phase 
is converged even at the lowest resolution considered here, 
128 3 . The linear growth phase of the instability terminates 
at t — 20, following which (v 2 ) 1 ^ 2 saturates and then de- 
cays for r > 30. During this phase of the evolution, (v 2 ) 1 ^ 2 
exhibits non-monotonic behaviour with resolution. As a re- 
sult, we conclude that {v 2 ) 1 ^ 2 is not a sufficient diagnostic 
to determine convergence in the non-linear regime. 

In the absence of explicit dissipative terms in the con- 
servation equations Q~H3] the effective (i.e., numerical) dissi- 
pation present in the simulation is governed by the choice 
of grid resolution. The numerical dissipation, expressed in 
units of diffusivity as (Ax) 2 / At, decreases with improved 
grid resolution for a fixed timestep. As grid resolution el- 
ements become finer, small-scale structures are preserved 



that would otherwise be smeared out by under-resolved sim- 
ulations whose numerical dissipation scale is too large to 
capture said structures. As small-scale structure can drive 
energy exchange and morphological evolution, particularly 
in the non-linear and turbulent regimes, establishing a con- 
verged solution is paramount for a physical interpretation 
of the simulation results. Expecting an exact point-to-point 
match of a quantity between different grid resolutions is in- 
appropriate given the turbulent nature of the problem at 
hand; therefore, we refer to convergence in the sense that 
quantities integrated over the entire volume do not change 
appreciably for a factor of two increase in N, the grid reso- 
lution in each dimension. 

Figure U shows the time evolution of the volume- 
averaged magnitude of the magnetic field for sets of 3D MHD 
simulations at various grid resolutions. This figure serves 
as a convergence study of the KHI simulations in the non- 
linear regime, as convergence is clearly obtained at N = 512. 
That is, the difference in evolution of magnetic energy in 
the non- linear regime is negligible between the N = 512 and 
N — 1024 simulations. Based on this figure and the discus- 
sion presented above, we conclude that the decay of turbu- 
lence in the non-linear regime is driven by physical, rather 
than numerical processes. We therefore treat = 512 as 
our fiducial resolution and the 3M512 run as our fiducial 
simulation. 

By contrast, 2D MHD simulations of the KHI do not 
exhibit convergence in the non-linear regime. Figure[S]shows 
the same quantity as in Figure [4] but for a series of 2D sim- 
ulations of the KHI at resolutions up to 16, 384 2 , with little 
evidence for convergence in the non-linear regime. This is 
evidenced by the absence of both a consistent peak in the 
magnetic field and a single sustained value at late times. 
The lack of convergence of 2D KHI simulations beyond the 
linear growth stage suggests that 2D studies of the non- 
linear and turbulent behaviour are inadequate; therefore, a 
full 3D treatment is required if one is interested in making 
quantitative statements beyond the linear growth stage. A 
more detailed comparison of the evolution of the two- and 
three-dimensional KHI is found in Figure H3 which shows the 
time evolution of the volume-averaged internal, kinetic, and 
magnetic energies, each normalised by the volume-averaged 
total energy. As can be seen from Figure |SJ the evolution of 
the energetics in the 2D system shows substantial differences 
from the 3D case. In particular, the 2D flow is more efficient 
than the 3D flow at generating magnetic energy from the 
available kinetic energy and this magnetic energy is less effi- 
ciently dissipated into heat. Internal energy is the dominant 
energy component and increases throughout the simulation 
because there is no cooling prescription. 

Further evidence of the differences in behaviour of the 
two- and three-dimensional KHI can be found through com- 
parison of Figures [Jj and [8] These figures show the time- 
averaged spectral distributions of the magnetic and kinetic 
energies in the three- and two-dimensional simulations, re- 
spectively, compensated by k A ^ to enable visual compari- 
son. Figure [7] shows that, for the 3D simulations, the spec- 
tral distribution of these quantities follow a fc~ 4//3 power-law 
for kL/(2n) > 3, over all the resolutions considered. The 
main effect of increasing resolution is to move the dissipa- 
tion scale to progressively smaller scales; from kL/(2ir) ~ 10 
(N = 128) to kL/(2-K) ~ 50 (N = 512). As evidenced by 
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Figure 7. Timc-avcraged, one-dimonsional spectral energy densi- 
ties for simulations 3M128 (dash-dot lines), 3M256 (dashed lines), 
and 3M512 (solid lines). Time averages are performed over the 
non-linear turbulent decay interval, [t pea k! *f]- Top panel: Kinetic 
energy power spectra, Ej^(k). Bottom panel: Magnetic energy 
power spectra, Eyi(k). Spectral energy densities are compensated 
by a factor of fc 4 / 3 . Consistent behaviour in both Epc(k) and 
Eyi(k) is seen in the non-linear turbulent decay regime across the 
resolutions studied, with an inertial range emerging at interme- 
diate scales for the 3M512 simulation. The effect of increasing 
resolution is to shift the dissipation scale to smaller scales (i.e., 
higher k). 
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Figure 8. Time-averaged, one-dimensional spectral energy den- 
sities for numerous 2D MHD simulations at different resolutions. 
Time averages are performed over the non-linear turbulent decay 
interval, [i pe aki *f]- Line colours are the same as in Figurc[5]and 
the resolution hierarchy can be deduced from the cut-off in k. 
Top panel: Kinetic energy power spectra, En(k). Bottom panel: 
Magnetic energy power spectra, Eyi(k). Spectral energy densities 
are compensated by a factor of A; 4 / 3 . With increased resolution, 
an approximately converged inertial range emerges in Em (k) and 
the dissipation scale moves to higher k. However, increasing res- 
olution alters the spectral distribution of kinetic energy, rather 
than merely pushing the dissipation scale to smaller scales. 



Figure [8] the behaviour of the 2D simulations is different. 
While the spectral distribution of the magnetic energy reser- 
voir shows approximate convergence to a k~ 4 ^ 3 power-law, 
the spectral distribution of the kinetic energy does not ap- 
pear to follow a simple power-law and shows a changing 
dependence as higher resolutions are probed. That is, the 
effect of increasing resolution in the 2D case is to alter the 
spectral distribution of kinetic energy, rather than simply 
move the dissipation scale to higher k as observed in the 3D 
case. This suggests the operation of an inverse cascad^3 m 
decaying 2D magnetised turbulence as has been known to 
occur in 2D hydrodynamic turbulence. In nature, turbulence 
must be inherently 3D and, as such, the clear differences be- 
tween the evolution of the KHI in two- and three-dimensions 
call into question conclusions reached regarding the appli- 
cability of the non-linear KH I to real physical systems from 
exclus i vely 2D studies ( e .g.. iFrank et all 1 19961; I Jones et al.l 
Il997bl ; Ijeong et al1l2000l : iBucciantini fc Del Zannall200rJ ). 



2 By "inverse cascade" we mean that energy in the magnetic 
reservoir is initially spectrally dominated at small scales and 
evolves to become primarily distributed on large scales. We do 
not mean to imply a dynamo process by using this phrase. 



5 EVOLUTION 

Comparing Figures [3] and [4] makes clear that significant 
magnetic field amplification only occurs within the non- 
linear stage of the instability (i.e., for r > 10). Figure [4] 
highlights that there are three regimes in the evolution of 
the magnetic field during the non-linear stage: amplifica- 
tion (10 < r < 30), saturation (30 < r < 50) and decay 
(t > 50). For purposes of clarity during subsequent dis- 
cussion, we further subdivide these regimes into intervals 
which are marked by vertical lines in Figure [4] and repre- 
sent, in chronological order, the times during the fiducial 
simulation, 3M512, at which the volume- averaged magnetic 
field magnitude relative to the initial field strength grows by 
10% (tw = 2.01), grows to |B pcak (tGi/3 = 3.15), grows to 
f-Bpeak (£g2/3 = 4.21), reaches B pca k (Vak = 6.43), decays 
to !-Bp ea k (£d2/3 = 9.86), and the time at the termination 
of the simulation (if = 15.00). The maximum amplitude 
attained by the magnetic field in the fiducial simulation is 
Speak = 7.27£>o an d the subscripts on the times are meant 
to indicate growth (G) and decay (D) stages of the magnetic 
energy. 

To visually assess the development of the KHI, slices of 
gas density, magnetic field magnitude, vorticity magnitude, 
and current density magnitude are shown in Figure for 
the fiducial simulation. At t = tio (not shown in Figure [9j|, 
the familiar linear growth wave-like pattern of the KHI is 
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Figure 9. 2D slices taken in the yz-plane at x = 0.5L from the 3M512 simulation at times *gi/3 (first row), t pe ^ (second row), and tf 
(third row). From left to right, the columns show the gas density, p, magnetic field strength relative to the initial value, B/ Bo, logarithm 
of the vorticity magnitude, log(|V X v|), and logarithm of the current density magnitude, log(| V X B|). The bracketed numbers above 
each column mark the [minimum, maximum] parameter values for the linear-scale color bar used to plot the respective quantity in all 
rows of the column. By the saturation time, t pea k> the initial shear layer is destroyed and the flow enters a state of decaying turbulence. 



developed, with numerous small-scale, low-pressure vortices 
forming that are the sites of magnetic field amplification 
due to twisting and stretching of field lines. Although ini- 
tially less pronounced than the small-scale vortices, a set of 
two large vortices along each shear layer begins to develop 
as a consequence of the single mode perturbation introduced 
into the computational box at to = 0. The KHI continues to 
develop into the non-linear stage by t = tGi/3, at which time 
two primary commensurate features have been established 
along each interface with multiple mini-vortices arising from 
secondary instabilities. The non-linear evolution continues 
and produces finger-like strands of density, magnetic field, 
and pressure by t = tG2/3 ( n °t shown in Figure |5J). When 
the magnetic field reaches its peak amplitude at t = i pe ak, 
the shear layer is nearly shredded beyond recognition into 



turbulence with evidence for the single mode form of the ini- 
tial perturbation also nearly erased. At this point in time, 
the magnetic energy production mechanism, namely, a con- 
tinually driven shear layer, is absent and the magnetic field 
begins to decay as the fluid motions remain turbulent and 
the fluid is well-mixed. Subsequent evolution of the system 
to late times is characterised by gradual decay of magnetic 
energy. 

The time-averaged kinetic and magnetic energy power 
spectra computed over these same time intervals from the 
fiducial simulation are shown in Figure I10I Table [3] pro- 
vides the spectral slopes for the intermediate scales, 5 < 
kL/(2n) < 30, corresponding to each time averaging inter- 
val in Figure I10I Figure [10] shows that magnetic energy is 
concentrated in small spatial scales as the KHI begins to 
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Figure 10. Time-averaged, one-dimensional spectral energy den- 
sities for the 3M512 simulation. From top to bottom are the power 
spectra for kinetic energy, magnetic energy, and ratio of magnetic- 
to-kinetic energies. Time averages are performed over the inter- 
vals [tio%i*Gl/3] {long-dashed lines), [tGi/3; *G2/3] (short-dashed 
lines), [t G 2/3:*pcak] [dash-dot lines), [f. pcak , t D2 /-s] (solid lines), 
and [tD2/3i*f] (dotted lines). As the simulation progresses, the 
spectral equipartition point shifts to larger scales, where mag- 
netic energy dominates kinetic energy on scales smaller than the 
equipartition scale. 

develop from the linear to saturated state (i.e., from t = t\o 
to t — tpcak). As the volume- averaged magnetic field is am- 
plified and peaks, magnetic energy at small scales (i.e., large 
k) grows with a fixed slope, while the spectral shape at 
small k flattens. By contrast, during this phase of the non- 
linear evolution, kinetic energy contained on large scales, 
kL/(2n) < 30, retains the same spectral slope and magni- 
tude, while the majority of the kinetic energy amplification 
occurs on small spatial scales, fci/(27r) > 30, due to the de- 
velopment of small-scale vortices. Once the peak magnetic 
energy is reached, magnetic energy on large scales decays 
more gradually than that on small scales, causing the spec- 
trum to steepen, while kinetic energy on large scales decays 
more rapidly than that on small scales, causing the spectrum 
to flatten. 

The bottom panel of Figure 1101 makes the comparison 
of the spectral evolution of magnetic and kinetic energies 
explicit. As the KHI develops from the linear regime to- 



E(k) 


m Gl/3 


Gl/3 
m G2/3 


G2/3 
peak 


peak 
m D2/3 


D2/3 
m. f ' 


E K (k) 
E M (k) 


-2.79(4) 
-0.31(2) 


-2.27(3) 
-0.44(3) 


-1.97(5) 
-0.80(2) 


-1.47(2) 
-1.33(2) 


-1.33(2) 
-1.62(2) 



Table 3. Slopes, m, of the kinetic, Ej^(k), and magnetic, E^(k), 
one-dimensional spectral energy densities from a log-log fit over 
the range, 5 < fcL/(27r) < 30, for the fiducial 3M512 simulation. 
Time averaging intervals for the spectral energy densities are de- 
noted by the subscript and superscript on m and conform to the 
notation described in the text (see Q. Uncertainties on the last 
significant digit are given in parentheses and correspond to the 
1-sigma level. At late times, both Ej^(k) and E]^(k) exhibit an 
inertial range approximated by a k~ 4 / 3 power-law. 

wards turbulence, the E M (k)/E-K(k) spectrum tends to in- 
crease and level off with increasing k. This behaviour was 
also observed in the relativ i stic M HD KHI simulations of 
IZhang. MacFadven fc Wanel (|2009l ). The equipartition point 
of magnetic and kinetic energies slides toward larger scales 
for the entirety of the evolution, until magnetic energy dom- 
inates over kinetic energy across nearly all scales at the ter- 
mination of the simulation. Although the individual kinetic 
and magnetic energy spectra are decreasing in amplitude 
after t = t pe ak, the equipartition point continues to shift 
towards lower k. 

The time-averaged transfer functions associated with 
energy exchange with the magnetic energy reservoir are 
shown in Figure [11] and provide a quantification of mag- 
netic energy sources/sinks as a function of scale k. Only the 
transfer functions associated with turbulent motions (i.e., 
the v t piece of the velocity decomposition) are plotted, as 
we found that the transfer functions associated with pure 
shear (i.e., 5*xyf) and cross terms (i.e., Xxyf) are negligible 
players in energy transfer in comparison. We first consider 
energy transfer during the stages of the KHI development 
leading up to saturation (i.e., from t = tio to t = £ pe ak)- The 
dominant growth mechanism of magnetic energy at large 
and intermediate scales is due to turbulent motions twist- 
ing/stretching magnetic field lines (i.e., TsKT(k) < and 
TkBT(fc) > 0). Transfer inside the magnetic energy reser- 
voir by turbulent velocities (i.e., TbbaCs)) is responsible for 
an inverse cascade of magnetic energy. Work done against 
magnetic pressure gradients by turbulent compressive mo- 
tions (i.e., Tbkp(&) and TkBp(fc)) is negligible in compar- 
ison to other magnetic transfer mechanisms. Although not 
plotted, we inspected the kinetic energy transfer functions 
and found the following behaviour. The dominant contribu- 
tion to large-scale kinetic energy growth between t = tio 
and t — tpeak is from advectioqj within the kinetic energy 
reservoir (i.e., Tkka(&) and Xkka(^)), with ancillary con- 
tributions from both compressible turbulent motions within 
the kinetic energy reservoir (i.e., TkKc(&) and XnKc(k)) 
and transfer from the internal energy reservoir by compres- 
sion (i.e., TjKc(k) and SiKc(k)). Meanwhile, energy is be- 
ing transferred out of the kinetic energy reservoir on these 
large scales into the magnetic energy reservoir by turbulent 

3 Here, and henceforth, "advection" is used to refer to the trans- 
fer of energy between scales but within the same form. For exam- 
ple, the "advection" of kinetic energy from large-to-small scales. 
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Figure 11. Time-averaged, one-dimensional spectral energy 
transfer functions associated with energy transfer into/out of the 
magnetic energy reservoir for the 3M512 simulation. From top to 
bottom are the transfer functions TbbaW> TbktWi TkbtWi 
^bkpW, and Tkbp(&)- Time averages are performed over the 
same intervals as in Figure 1101 with the same line style conven- 
tion used. The exact details of the energetics are seen to be 
highly time-dependent. In general, the transfer function ampli- 
tudes evolve in time, but their spectral shape remains fairly con- 
sistent, albeit with some translation in k. 
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Figure 12. Time-averaged, one-dimensional spectral energy 
transfer functions associated with energy transfer into/out of 
the magnetic energy reservoir for the 3M512 simulation for the 
time-averaging interval [f pea k;if]- The fluid is in a state of de- 
caying turbulence over this interval in time. The lines shown 
here are a representation of the data from the solid lines and 
dotted lines in Figure 1111 but placed on the same scale to al- 
low for relative comparisons. The transfer functions shown are 
Tbba(^) (solid line), TbktCs) (dotted line), Tkbt(&)i (dash-dot 
line), TbkpWi (short-dashed line), and Tkbp(&)i (long-dashed 
line) . Energy transfer is dominated by exchanges within the mag- 
netic energy reservoir (i.e., TbbaCO) an d transfer mediated by 
magnetic tension (i.e., Tbkt(^) and Tkbt(^))- Magnetic pres- 
sure effects are small in comparison due to the weakly compressive 
nature of the subsonic and sub-Alfcnic flow studied here. 



fluid motions twisting/stretching the magnetic field (i.e., 
TBKT(fc) < 0). The small-scale growth of kinetic energy dur- 
ing this time is overwhelmingly dominated by the same mag- 
netic tension forc<|3 acting on turbulent motions that causes 
kinetic energy loss on large scales. In summary of the KHI 
evolution leading up to saturation, we find that the magnetic 
field grows first at small scales and then cascades to larger 
scales, which is evidence for an inverse cascade operating 
in the KHI, and the dominant energy exchange mechanism 
involves turbulent fluid motions interacting with magnetic 
tension. 

We now turn our attention to times after t — £ pe ak, 
where the fluid is in a turbulent state and energy decays 
away by numerical dissipation. From t — t poa k onward, when 
the simulation box is fully turbulent and the shear layer is 
destroyed, an abrupt transition occurs where the subsequent 
evolution of the kinetic energy spectrum over all scales is 
determined primarily by interactions with the magnetic en- 
ergy reservoir. Figure [12] shows the energy transfer between 



4 Here, and henceforth, "tension" is used to describe the restoring 
force directed along the radius of curvature that is exerted by bent 
magnetic field lines. We do not mean to imply that the magnetic 
field is always putting tension on the fluid. 
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Figure 13. Diagram showing the contributions from transfer functions to the two-way energy exchange between the kinetic and magnetic 
energy reservoirs at late times. Percentages in parentheses indicate the amount of energy exchange into that reservoir as described by 
the associated transfer functions, relative to the time-averaged total energy exchange rate, (tJjSJ?'), over the time interval [t pea k> if]- 
The dominant scales (i.e.. small, large, all) across which the energy transfer operates are indicated for each exchange path. For instance, 
transfer of large-scale kinetic energy into the magnetic energy reservoir by twisting/stretching of magnetic field by fluid motions (i.e., 
Tbkt < and Shut < 0) is responsible for 40% of the total energy exchange. The line styles are chosen to overlap with those of Figure 
1121 Magnetic tension is the dominant transfer mechanism for exchanges into/out of the magnetic energy reservoir. The kinetic/magnetic 
energy reservoir interactions result in a net magnetic energy gain rate. This energy then cascades from large-to-small scales and is further 
exchanged forwards and backwards with the kinetic energy reservoir until it is ultimately dissipated. 



magnetic energy transfer functions in the time-averaged de- 
cay stage from t = t pca k to t = tf. As before, exchanges 
within the magnetic energy reservoir and transfer mediated 
by magnetic tension dominate the magnetic energetics, with 
exchange by magnetic pressure gradients being negligible 
across all scales. Magnetic energy is supplied by large-scale, 
kL/(2n) < 10, kinetic energy loss due to turbulent fluid 
motions working against the magnetic tension force, as ev- 
idenced by negative values of Tbkt (k) on large scales. The 
positive values of Tbkt(&) on intermediate scales peaking 
at kL/ (2tv) ~ 80 indicate that magnetic energy is also being 
placed into intermediate-scale kinetic energy by the reversal 
of the process just described. Note that the transfer func- 
tion Tbkt(&) reveals that a significant amount of energy is 
being exchanged, but one cannot say on what scales it is dis- 
tributed in the magnetic energy reservoir. Presumably, some 
of this large-scale kinetic energy is transferred into large- 
scale magnetic energy. The kinetic energy reservoir con- 
tributes a modest amount of intermediate/small-scale mag- 
netic energy via work done on the magnetic field by fluid 
motions, as shown by positive values of Tkbt(&)- Most of 
the small-scale magnetic energy comes from turbulent trans- 
fer within the magnetic energy reservoir from large-to-small 
scales (i.e., Tbba transitions from negative to positive val- 
ues going from large-to-small scales). Thus, Figure [P2l tells 
a story of a mechanism for ongoing large-scale magnetic en- 
ergy production and a turbulent cascade from large-to-small 
scales within the magnetic energy reservoir. This small-scale 
energy is exchanged forwards and backwards with the ki- 
netic energy reservoir and is gradually dissipated, allowing 
the magnetic field to keep a relatively sustained value in the 
absence of a driven shear layer. 

We note that in the ideal MHD limit (i.e., in the ab- 
sence of numerical dissipation) there is no mechanism for 



Transfer Rate T B kt Sbkt Tkbt Skbt Tbkp Sbkp Tkbp 
(%) (%) (%) (%) (%) (%) (%) 

{dE+/dt) 53.9 7.4 29.0 6.8 0.5 2.4 0.0 

(dE^/dt) 69.5 0.6 11.7 0.0 5.0 4.6 8.6 



Table 4. Percentage breakdowns of the contributions to the mag- 
netic energy gain and loss rates by transfer function. The leftmost 
column is the time-averaged magnetic energy transfer rate, where 
the superscripts denote + for gain and — for loss. The right- 
ward columns list the percentage contribution from each transfer 
function involved in magnetic energy exchange. Notably, energy 
transfer mediated by turbulent motions (i.e., the vt component 
of v) interacting with the magnetic tension force are the primary 
players in energy exchange with the magnetic energy reservoir. 

the direct transfer of internal energy to magnetic energy be- 
cause there is no internal energy contribution in the induc- 
tion equation; therefore, the magnetic and internal energy 
reservoirs cannot engage in energy exchange. Perhaps un- 
surprisingly, the evolution of the internal energy spectrum is 
determined by transfer between the kinetic and internal en- 
ergy reservoirs by advective and compressive motions, while 
numerical dissipation effects will also change the internal 
energy. 

Finally, we summarise the nature of energy transfer 
operating in the KHI over the late-time turbulent decay 
stage from t pea k to tf. Separately collecting the positive 
and negative contributions from each transfer function al- 
lows one to determine the total magnetic energy gain and 
loss rates via exchanges described by the transfer func- 
tions. We find a time-averaged magnetic energy gain rate 
of (dE^/dt) = 4.4 x and loss rate of (dE^/dt) = 

—2.3 x 10~ 3 , which gives a net magnetic energy gain rate 
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due to energy transfer with the magnetic energy reservoir 
of (dE^f /dt) = 2.1 x 10~ 3 , all in code units. Note that 
the transfer functions describing the magnetic energy cas- 
cade (i.e., Tbba(&) and iSbba(&)) are not included in this 
inventory because they cannot contribute to overall mag- 
netic energy gain or loss. Simultaneous with the magnetic 
energy reservoir exchange described by the transfer func- 
tions are net magnetic energy loss rates resulting from both 
the decaying nature of the MHD turbulence and numeri- 
cal effects. The time-averaged magnetic energy decay rate is 
(dEm/dt) — —1.6 x 1CP 3 and the time-averaged loss term 
associated with numerical dissipation of magnetic energy is 
(Dm) = -3.0 x 1CT 3 . Table g] lists the relative contribu- 
tions of each transfer function involved in magnetic energy 
gain and loss rates averaged from £ poa k to tf. Stretching and 
twisting of magnetic field lines by the turbulent velocity 
field (i.e., Tbkt(&) and TkBT(fc)) is the dominant exchange 
mechanism at work during late times, accounting for 83% 
and 81% of energy transfer leading to magnetic energy gain 
and loss, respectively, while magnetic pressure is a negligi- 
ble contributing transfer mechanism for the subsonic and 
sub-Alfenic flows we consider. 

To understand the two-way energy flow into/out of 
the magnetic energy reservoir, we consider the total time- 
averaged magnetic energy exchange rate, (dEffi/dt) = 



UdE+/dt) 



+ (dE M /dt) , and construct a schematic diagram 



in Figure [13] that tracks the contributions of each transfer 
function to (dEffi/dt). Figure \T3\ illustrates that the kinetic 
energy reservoir interacts with the large-scale field and in- 
jects energy into the magnetic energy reservoir. This energy 
cascades down to smaller scales and is exchanged backwards 
and forwards with the kinetic energy reservoir, before ulti- 
mately being dissipated. The turbulent cascade from large- 
to-small scales (i.e., Tbba(&)) operates on 61% of (dE^/dt), 
making the cascade within the magnetic energy reservoir an 
effective mechanism for breaking down magnetic structures. 



6 DISSIPATION 

The extremely large Reynolds numbers that characterise as- 
trophysical flows suggest that it is appropriate to carry out 
numerical simulations of these same flows in the inviscid 
flux-freezing regime, where explicit dissipation terms are 
omitted from the momentum and induction equations. In 
nature, however, astrophysical flows have some small, but 
finite amount of viscosity and resistivity, which violates the 
assumption of an inviscid, flux-frozen flow. In a turbulent 
flow, such as that considered here, it is these dissipation 
terms that mediate the dissipation of small-scale turbulent 
structures and conversion of magnetic and kinetic energy 
contained in these structures into thermal energy. When 
performing calculations in the inviscid, flux-freezing regime, 
simulators hope that the details of dissipation, which are 
provided by the algorithm, have little influence on large-to- 
intermediate scales. If dissipation does influence these scales, 
simulators hope that the details of the numerical dissipation 
are sufficiently similar to physical dissipation such that the 
simulation remains an accurate representation of the phys- 
ical system. This non-trivial issue regarding the validity of 
relying on numerical dissipation to adequately capture the 
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Figure 14. Timc-averagcd, onc-dimcnsional, shell-averaged spec- 
tral energy densities for simulations with (red lines) and without 
(black lines) explicit dissipation. Simulations with explicit dis- 
sipation introduced at saturation are: 3M512D^ (long-dashed 

lines); 3M512D^ (solid lines); 3M512D^^ (dash-dot lines); 
3M512D^ (triple dot-dash lines); 3M512Dj^ (dotted lines). 
Simulations without explicit dissipation are: 3M512 (solid lines); 
3M256 (dashed lines); 3M128 (dotted lines). Time averages are 
performed over [t pe ak; tf]. Top panel: Kinetic energy power spec- 
tra, En((k)). Bottom panel: Magnetic energy power spectra, 
Eyi((k)). Both spectral energy distributions are compensated 
by fc 4 / 3 . For all dissipation coefficients explored, power is de- 
pleted on small scales for the runs with explicit dissipation com- 
pared to the fiducial ideal MHD run. Simulations 3M512^g^ and 
3M512j^, which incorporate explicit dissipation terms, provide 
close matches to the ideal MHD simulations 3M512 and 3M256, 
respectively. 



behaviour of physical dissipation in simulations is what we 
address in the present section. 

Figure [14] compares shell-averaged (see fj3j spectral en- 
ergy densities obtained from simulations of decaying turbu- 
lence arising from the KHI with explicit dissipation added 
to the momentum and induction equation (specifically, we 
include the effects of kinematic shear viscosity and Ohmic 
resistivity) compared to the data of the convergence study 
presented in Figure [7] The simulations with explicit dissi- 
pation were initialised from the fiducial simulation, 3M512, 
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Figure 15. Time-averaged, one-dimensional, shell-averaged spec- 
tral energy densities for simulations with (red lines) and without 
(black lines) explicit dissipation. Simulations with explicit dis- 
sipation introduced at saturation are: 3M512D^ (solid lines); 
3M256D^ (dashed lines); 3M128D^ (dash-dot lines). Simula- 
tions without explicit dissipation arc: 3M512 (solid lines); 3M256 
(dashed lines); 3M128 (dash-dot lines). Time averages are per- 
formed over [tpeak> if]- Top panel: Kinetic energy power spectra, 
Eyi((k)). Bottom panel: Magnetic energy power spectra, £Jjy[((fc)). 
Both spectral energy distributions are compensated by fc 4 / 3 . Nu- 
merical dissipation becomes a more dominant contributor to the 
total dissipation with decreasing numerical resolution. 

at t = tpcak- We found that kinematic shear viscosity and 
Ohmic resistivity coefficients, v = 3.25 x 10~ 6 and r\ — 
2.125 x 10 -6 (Pm = v/rj = 1.53, simulation 3M512D^), 
produced a small, but non- negligible change in the mag- 
netic and kinetic spectral energy densities compared to the 
fiducial simulation. The dissipation coefficients were then 
increased by factors of two (at fixed magnetic Prandtl num- 
ber, Pm, and initialised from t — £ pea k of 3M512) until the 
magnetic and kinetic spectral energy densities provided a 
close match to those obtained from simulation 3M256. This 
occurs for simulation 3M512D^, where = 2.6 x 10 -5 
and ?7fid = 1-7 x 10~ 5 (where subscript 'fid' denotes that we 
treat these as our fiducial values), a factor of 8 increase over 
the dissipation coefficients found to match 3M512. This re- 
sult suggests that the numerical dissipation present in the 
ideal simulations scales as (Aa;) 3 rather than (Aa;) 2 as would 



be expected for second-order convergence. While the algo- 
rithms in Athena are overall second-order, the spatial re- 
construction method used in the ideal simulations is third- 
order, perhaps suggesting that the numerical dissipation in 
the KHI problem is determined by the order of spatial recon- 
struction. To demonstrate that this scaling holds generally 
for Athena-run KHI models, however, would require an en- 
semble of 3D simulations including dissipation and, which is 
beyond the scope of this work. 

Figure [TS] examines the convergence of the simulations 
using the fiducial dissipation coefficients, where the data 
from the ideal MHD convergence study presented in Fig- 
ure[7]are included for comparison. At resolutions lower than 
N = 512 (i.e., N = 128, 256) we find that numerical dissipa- 
tion plays an increasingly important role. In particular, there 
is a close correspondence across all scales between the simu- 
lations at N = 128 with (3M128D^) and without (3M128) 
contributions from explicit dissipation, indicating that so- 
lutions at this resolution are dominated by effects due to 
numerical dissipation. The TV = 256 case with the fiducial 
dissipation coefficients, 3M256D;[^, matches the large-scale 
behaviour of 3M128 and the small-scale behaviour of 3M256. 
As noted previously, the N = 512 case with the fiducial dissi- 
pation coefficients, 3M512DJJJ, provides a close match to re- 
sults obtained for 3M256 at all scales. The primary difference 
from 3M256 is a small power deficit for 3M512D^ at scales 
around the dissipation scak0, kL/(2ir) ~ 30, due to power 
being transferred over to smaller scales, kL/(2n) ~ 100, 
where there is a slight power excess. We regard the sim- 
ulation using the fiducial dissipation coefficients as being 
converged at N = 512 for two reasons. First, we already 
demonstrated that simulations conducted in ideal MHD are 
converged at this resolution (see Q, implying that numeri- 
cal dissipation plays a small role in simulations at this reso- 
lution. Second, for simulations incorporating physical dissi- 
pation, convergence implies that the dissipation scale is re- 
solved. As elucidated above, the location of the dissipation 
scale for simulation 3M512D^ has moved to smaller k (i.e., 
larger physical scales) than is the case in the ideal simula- 
tion at this same numerical resolution, 3M512. This implies 
that the location of the dissipation scale is determined by 
the dissipation terms themselves rather than numerical ef- 
fects and as a result, we can conclude that the dissipation 
scale associated with i/Rd, ??fld is resolved at N = 512. 

With these arguments in mind, Figure [16] compares 
transfer functions associated with energy exchange with the 
magnetic energy reservoir for simulations 3M512D^ and 
3M256, time-averaged over the interval [i pe ak,if]- At large 
spatial scales, k < 5, the transfer functions for 3M512DJ5J 
and 3M256 are well-matched. At intermediate scales, 5 < 
k < 30, the transfer functions are well-matched for transfer 
from magnetic energy to magnetic energy through advec- 
tion, Tbba(&), and from magnetic energy to kinetic energy 
through tension forces, Tbkt(&). By contrast, we see greater 
transfer from kinetic to magnetic energy through tension, 
TkBT(fc), at these scales for the 3M512D^ simulation. At 



5 In this work, the dissipation scale refers to the approximate 
turn-over scale where the spectral energy distribution transitions 
from a power-law inertial range to the steep decline at small 
scales. 
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Figure 16. Time-avcragcd, onc-dimcnsional spectral energy 
transfer functions associated with energy transfer into/out of the 
magnetic energy reservoir for the ideal MHD run 3M256 (black 
lines) and the simulation 3M512D]^ (red lines) where explicit 
dissipation was introduced at saturation with fiducial dissipation 
coefficients. The time-averaging interval and choice of line styles 
are consistent with those of Figure 1121 The dominant transfer 
functions in energy exchange (Tbba(^) an d TbktC 2 )) are g en_ 
erally well-matched between these two runs at scales larger than 
the dissipation scale (i.e., k < 30), suggesting that the physics 
of energy transfer in MHD turbulence is robust to the effects of 
numerical dissipation. 



small scales, k > 30, we see that peaks in the transfer func- 
tions are shifted to smaller scales in 3M512D1JJ , as compared 
to 3M256, as a consequence of the higher numerical resolu- 
tion in this simulation. Finally, at all scales, the effect of 
dissipation is to reduce the (already small) contribution of 
energy transfer through compressive motions, Tkbp(&) and 
7bkp(&). Overall, these results demonstrate the robustness 
of the physics of energy transfer within decaying MHD tur- 
bulence to the effects of numerical dissipation at scales larger 
than the dissipation scale. 

The ideal MHD equations are modified in the presence 
of explicit dissipation by including the viscous term, i/(V-t), 
in the momentum equation and the resistive term, ?7V 2 B, 
in the induction equation, where v and r\ are the assumed 
constant kinematic viscosity and Ohmic dissipation, respec- 
tively, and r is the stress tensor for an isotropic fluid. The 
transfer functions associated with explicit dissipation take 
the form (see Appendix |B|) . 



T„{k) = [pv] 



- (V • t) 
P 



(k) + v- (V 



(*) 



T n {k) 



B 



V 2 B (k) 



(21) 
(22) 



These are plotted in Figure [17] for simulations 3M256 and 
3M512D}^. Note that for ideal (i.e., inviscid and non- 
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Figure 17. Time- averaged, one-dimensional dissipation trans- 
fer functions for the 3M256 (black lines) and 3M512D}^ (red 
lines) simulations over the interval, [t pea k! *f] an d compensated 
by — fc -2 / 3 . Top panel: Kinetic energy dissipation transfer func- 
tions and numerical dissipation terms. Bottom panel: Magnetic 
energy dissipation transfer functions and numerical dissipation 
terms. See the legends within each panel for the specific quanti- 
ties being shown. 



resistive) MHD simulations, such as 3M256, the physical 
dissipation transfer functions T v (k), T v (k) of course do not 
contribute to the overall energy transfer inventory, but are 
instead computed for the sake of establishing 'effective' 
quantities. The effective dissipation transfer function data 
of 3M256 adopt ^fid,?7fid to enable comparison with T v (k), 
T n (k) from 3M512D^. Also shown in Figure [T71 for both 
simulations are the quantities _Dk(&), -Dm(&), derived by 
calculating the residual of the terms in Equations 1 1 71 and 1181 
(including T v {k) and T v (k) for 3M512D^) and the total dis- 
sipation (£ K (fc) = T v (k) + D K {k), f M (fc) = T v (k) + D M {k)) 
for simulation 3M512D 1 ^. Figure [171 shows that the spec- 
tral distribution of T v (k), T v (k) is very similar between the 
two simulations for scales larger than the dissipation scale 
and that the spectral distribution of numerical dissipation 
in 3M256 is very close to that expected from physical dissi- 
pation close to the grid scale (i.e., small scales). A further 
point comes from comparing physical and numerical dissi- 
pation in simulation 3M512D]^. For this simulation, we find 
that, while physical dissipation, T v (k), dominates over nu- 
merical dissipation, Dyi{k), for the magnetic energy dissi- 
pation terms by a factor ~ 100 for 2 < kL/(2n) < 30, the 
same is not true for the kinetic energy dissipation terms, 
where physical dissipation, T v (k), and numerical dissipa- 
tion, Dj((k), are relatively comparable to within a factor 
of ~ 2 over this range of scales. Similar levels of numerical 
and physical kinetic energy dissipation could be due to the 
computation of derivatives that are required for the viscous 
stress tensor, t, and the associated divergence, V • r. The 
same considerations do not apply for the addition of Ohmic 
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Figure 18. Time-averaged, one-dimensional spectral energy den- 
sities for simulations where physical dissipation was incorpo- 
rated from the very beginning of the KHI simulation. Time av- 
erages are performed over the late-time non-linear turbulent de- 
cay interval, [i pGa k, if]- Top panel: Kinetic energy power spectra, 
Ej^(k). Bottom panel: Magnetic energy power spectra, Ejy[(fc). 
Both spectral energy distributions are compensated by fc 4 / 3 . The 
results from simulations with dissipation coefficients (y, n) = 
(2, 2), (1, 1), (2, 1), (1, 2) are shown by the dashed, solid, dotted, 
and dash-dot red lines, respectively. For reference, the ideal MHD 
simulations 3M512 (fiducial) and 3M256 are shown by the solid 
and dashed black lines, respectively. Enforcing explicit dissipation 
at the start of the KHI simulations tends to transfer small-scale 
power to large/intermediate scales as compared to the spectral 
distributions from ideal MHD simulations; however, no simple 
pattern emerges from changing the magnetic Prandtl number. 



diffusion to the induction equation due to the use of the CT 
algorithm for these terms, which may explain why physical 
magnetic energy dissipation greatly exceeds numerical mag- 
netic energy dissipation on scales larger than the dissipation 
scale. 

Finally, we consider the impact of changing the mag- 
netic Prandtl number on the physics of decaying turbulence. 
Simulations conducted in ideal MHD can consider one Pm, 
which is intrinsic to the algorithm. We measured Pm — 1.5 
for the application considered here. In nature, however, mag- 
netised flows can occur at a range of Pm; therefore, it is 
important to assess the impact that varying the magnetic 
Prandtl number can have on the flow. Figure [18] shows the 
change in the magnetic and kinetic spectral energy distribu- 
tions for 0.75 < Pm < 3.1, where these runs incorporated 
explicit dissipation from the very start of the KHI simulation 
(i.e., t = to). The apparent effect of changing the magnetic 
Prandtl number is to transfer power from small to interme- 
diate/large scales within the turbulence; however, any such 
trend does not appear to follow a simple linear pattern and 
so a clear measurement would require extensive further sim- 
ulation, which is beyond the scope of this work. 



7 SUMMARY AND DISCUSSION 

We performed a suite of 2D and 3D simulations of the KHI 
in the weakly magnetised, subsonic regime with a non-driven 
shear layer, focusing on the results of a high-resolution 3D 
MHD simulation. The problem setup, though simple and 
straightforward, was scrutinised in detail, paying particu- 
lar attention to dimensionality (2D vs. 3D), convergence, 
and properly resolving the shear layer in order to make a 
convincing argument for the physical nature of the KHI de- 
velopment beyond the linear growth. After establishing the 
basic evolution of energetics using volume-averaged energies 
and time-averaged energy power spectra, we took advantage 
of the energy conserving nature of Athena to investigate the 
spectral structure of the KHI development into MHD turbu- 
lence using spectral energy transfer function analysis. This 
analysis was then extended to characterise both numerical 
and physical dissipation in Athena. Here, we discuss our re- 
sults. 

Two-dimensional MHD simulations of the KHI (e.g. . 



Frank et all Il99tj; IJones et all Il997bl ; IJeong et all l200d ; 



Bucciantini fc Del Zannal 120061 ) are attractive due to their 
ability to achieve high resolution relative to their 3D coun- 
terparts; however, a demonstration of convergence of the 
resulting turbulent flow is required to justify 2D studies of 
MHD turbulence arising from the KHI. We observe well- 
converged solutions of the initial growth of the 2D KHI 
at the moderate resolution 512 2 , which justifies the linear 
growth stage of the 2D KHI as a hi ghly reliable, robust test 
for co de verification as suggested bv lMcNallv. Lyra fc Passvl 
(2012). However, the saturated state and level of magnetic 
energy sustainment fails to converge even out to the ex- 
tremely high-resolution 16, 384 2 , as evidenced by the time 
evolution of the volume-averaged magnetic field strength 
(see Figure[S]) and the changing shape of spectral energy den- 
sities with resolution (see Figure [5]). In stark contrast to the 
2D case, 3D simulations of the KHI are reliably converged at 
a resolution 512 3 over the full-course of evolution out to the 
turbulent and decaying stages. Furthermore, any real astro- 
physical scenario involving the magnetised non-linear KHI 
will certainly develop 3D character. Therefore, not only will 
2D treatments fail to capture the relevant physics due to 
convergence issues beyond the linear growth stage, they also 
misrepresent the inherently 3D nature of the astrophysical 
situations to which they are being applied. 

Time evolutions of volume-averaged energetics and 
slices of the simulation volume reveal the growth of magnetic 
energy to a saturated level at which time the shear layers 
are almost completely disrupted. The subsequent evolution 
leads to turbulence with a sustained, but gradually decay- 
ing, magnetic field. This general evolution is also observed 
in relativistic MHD simulati ons of the KHI when the driving 
mechanism is switched off (IBucciantini fc Del Zannal 120061; 

120091 ; IZrake fc MacFadvenl 
either a discontinuous shear 



Zhang. MacFadven fc Wand 
20 111 ). These studies adopt 
layer, use a Riemann solver of type HLLE, or both, yet 
we find that the sharp decline in kinetic energy during the 
non-linear growth and generation of a sustained magnetic 
field is robust to the details of the initial setup and Riemann 
solver used (see Appendix [A")) . We confir m the results of the 
relativ istic MHD study of the KHI of iBeckwith fc Stonel 
(|201ll ) in the Newtonian regime using a linearised Riemann 
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Figure 19. Ratio of total (i.e., numerical + physical, if applica- 
ble) magnetic-to-kinetic energy dissipation rates for simulations 
3M512 (solid line), 3M256 (dashed line), and 3M512D^ (dotted 
line), time-averaged over the interval [t pea k! tf]. The horizontal 
solid line marks the cquipartition of total magnetic and kinetic en- 
ergy dissipation rates. Magnetic energy dissipation exceeds kinetic 
energy dissipation for intermediate-to-small scales (i.e., k > 10). 



solver. While the generic result of the appearance of a 
saturated state is unaffected, we caution against using 
a setup with an unresolved interface and/or the HLLE 
Riemann solver for quantitative studies of the KHI. 

The spectral distributions of kinetic and magnetic en- 
ergy for 3D KHI simulations at late times follow an ap- 
proximate fc -4 ' 3 power-law on intermediate scales, 5 < 
kL/(2n) < 30, remaining unaltered for all resolutions con- 
sidered (see Figure [T}. A spectral slope oc k~ 4 ^ 3 over in- 
termediate scales also appeared in the strong-field driven 
supers onic MHD turbulence studies of iLemaster fe Stone] 
(|2009T > for resolutions of 512 3 and 1024 3 . The effect of in- 
creasing numerical resolution is to move the dissipation scale 
to smaller scales. The magnetic-to-kinetic energy spectral 
equipartition point shifts to larger scales throughout the 
simulation evolution (see Figure 1 10p . Performing a study 
of relativistic, ideal MHD turbulen ce arising from the KHI, 
IZhang, MacFadven fc Wand (120091 ) claim that this observed 
evolution of the Eyi{k) / E^(k) equipartition point indicates 
that the kinematic viscous dissipation is more efficient than 
the magnetic resistive dissipation; however, no direct study 
of dissipation was performed. Figure \T§\ shows the ratio of 
total magnetic-to-kinetic energy dissipation rates in the tur- 
bulent regime for simulations with (3M512D}^) and with- 
out (3M256, 3M512) explicit dissipation included. Figure [19] 
demonstrates that magnetic energy dissipation actually ex- 
ceeds kinetic energy dissipation across the majority of scales, 
k > 10. Therefore, the shift in the Eyi{k) / E^Qi) equiparti- 
tion point in Figure [10] is instead a consequence of the ex- 
change of large-scale kinetic energy into the magnetic energy 
reservoir mediated by turbulent motions acting against mag- 
netic tension (i.e., fluid motions twisting/stretching mag- 
netic field lines), as evidenced by the dominating negative 
values the TBKT(fc) in Figures llll and 1 1 21 The ambiguity 
as to whether dissipation rates or large-scale kinetic energy 
loss to the magnetic energy reservoir is the true mechanism 
behind the shift in the EM(k)/Ejc(k) equipartition point to 
large scales was resolved by appealing to transfer function 
analysis, illustrating that this is a powerful tool for studying 
how energy is transferred across scales and forms. 

Spectral energy transfer analysis allows for both the 
scale-by-scale quantification of energy transfer between 



reservoirs and identification of the mechanism responsible 
for the energy exchange, which are inaccessible from power 
spectra alone. As the KHI develops to a saturated state, the 
growth of magnetic energy is dominated by the magnetic 
tension force interacting with turbulent motions and an in- 
verse cascade is observed, meaning that magnetic energy is 
initially concentrated on small scales and then evolves to 
a spectrum dominated on large scales (see Figure 11 1 pi . At 
late times following saturation when the fluid is in a decay- 
ing turbulent state, we find no evidence for dynamo oper- 
ation for a single fluid treatment, as claimed from simula- 
tions of dec aying turbulence arising from relat ivistic MHD 
KHI studies lzhang, MacFadven fc Wans! (|200gh . Kinetic en- 
ergy contained in turbulent fluid motions is transferred to 
magnetic energy, primarily mediated by interactions with 
the magnetic tension force, and a turbulent cascade from 
large-to-small scales operates within the magnetic energy 
reservoir. This small-scale magnetic energy is interchanged 
forwards and backwards with the kinetic energy reservoir 
and is eventually dissipated, allowing the magnetic energy 
to decay. For the subsonic and sub-Alfvenic relative flow 
considered in this work, compressible effects are of ancillary 
importance in energy transfer. 



By their nature, numerical simulations exhibit dissipa- 
tive behaviour due to finite numerical resolution. Even in 
instances where physical dissipation terms are explicitly in- 
cluded in the solution of the MHD conservation equations, 
numerical dissipation is still present at some level. We found 
that the most important effect of increasing numerical res- 
olution for ideal MHD simulations was to move the dissi- 
pation scale to progressively smaller scales. While energy 
dissipation in ideal MHD simulations occurs preferentially 
on the grid scale, physical dissipation should act across all 
scales; therefore, determining the extent to which numerical 
dissipation affects MHD turbulence when physical dissipa- 
tion is present is nontrivial. We found that when the numer- 
ical resolution was held fixed, the location of the dissipation 
scale moves to larger spatial scales when physical dissipation 
is incorporated (3M512D 1 ^) compared to the corresponding 
ideal MHD simulation (3M512). This result indicates that it 
is the dissipation terms that determine the dissipation scale, 
rather than numerical effects. The physical dissipation scale 
is considered to be resolved when the dissipation scale (i.e., 
the turnover in the power spectrum at large k) moves to 
larger spatial scales than in the case without explicit dissi- 
pation terms included. In this sense, the "effective" resolu- 
tion of the simulation, by which we mean the location of the 
dissipation scale, is reduced by construction. Furthermore, 
when physical dissipation is introduced, the magnitude of 
numerical dissipation is diminished and the spectral charac- 
ter of the transfer functions (i.e., general shapes and relative 
proportions) involved in exchange with the magnetic energy 
reservoir are well-matched to their ideal MHD counterparts. 
These observations indicate the robustness of the physics of 
energy transfer in decaying MHD turbulence to the effects 
of numerical dissipation, at least for scales larger than the 
dissipation scale where numerical effects do not dominate. 
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8 CONCLUSIONS 

In light of the summary and discussion given in 5j7] our con- 
clusions are listed here as a bulleted list followed by some 
astrophysical implications of this work. 

• 2D KHI simulations do not converge beyond the linear 
growth stage, suggesting that application of these calcu- 
lations to astrophysical objects is inappropriate. 

• 3D KHI simulations reliably converge across all stages of 
evolution. The main effect of increasing numerical resolu- 
tion is to push the numerical dissipation scale to smaller 
spatial scales without changing the shape of the power 
spectrum. 

• For subsonic, weakly magnetised, decaying turbulence 
arising from the KHI, the spectral distributions of kinetic 
and magnetic energy for 3D simulations follow an approx- 
imate fc~ 4//3 power-law on intermediate scales. 

• We find no evidence for the operation of a dynamo in a 
single fluid system. At late times corresponding to decay- 
ing MHD turbulence, energy is injected into the magnetic 
reservoir as a result of kinetic energy interactions with 
the large-scale magnetic field. This magnetic energy tur- 
bulently cascades down to smaller scales and is exchanged 
backwards and forwards with the kinetic energy reservoir, 
before ultimately being dissipated. 

• The addition of physical dissipation has little impact on 
the overall shape of the kinetic and magnetic energy power 
spectra, even when the dissipation scale is well-resolved. 

• Incorporating explicit dissipation terms reduces the im- 
portance of numerical dissipation and moves the dissipa- 
tion scale to larger spatial scales. 

• The nature of numerical dissipation does not affect the 
physics of energy transfer within decaying MHD turbu- 
lence at scales larger than the dissipation scale. 

Our investigation of the subsonic KHI in the weak mag- 
netic field limit and the generalised spectral energy transfer 
function techniques we exploit serve as a launching point 
for future studies of MHD turbulence and the extension to 
more targeted astrophysical applications of the KHI. 

Non-ideal MHD effects have been observed to dramati- 
cally alter the non-linear behaviour of the KHI i n 2D multi- 
fluid MHD simulations (| Jones fc Downed lioill ) . Most no- 
tably, the Hall effect, when dominant, is seen to cause a 
dynamo effect that leads to a perpetually growing magnetic 
field in the non-linear regime, meaning that the KHI never 
achieves a truly saturated state when n on-ideal MHD effects 
are considered (|jones fc Downesll201ll ). This could have sig- 
nificant relevance to the physics of the KHI in weakly ionised 
astrophysical plasmas, as expected in protoplanetary accre- 
tion discs and molecular clouds. We demonstrated that even 
ideal MHD simulations of the 2D KHI do not converge to 
a saturated state; therefore, 3D non-ideal MHD simulations 
are required to test the suggestion that the Hall effect leads 
to a dynamo resulting in non-convergence of the saturated 
state. 

In addition to serving as a direct examination of 
KHI physics, this work provides a valuable baseline for 
investigations of shear layers in astrophysical systems 
also subject to the family of current-driven instabilities 
(CDI). Such systems could potentially fea ture either sharp 
shear layers, such as those explored by iBatv fc Keppensl 



(|2002h or iMizuno. Hardee fc Nishikawal (|201ll ). or more 
gradual profiles, such a s tho se examined analytically by 
iNalewaiko fc Begelmanl i|2012T l. Regardless of the details, 
it should be possible to compare growth rates of systems 
unstable to the KHI and CDI to both analytic estimates 
of linear growth and those rates measured empirically in 
this work. This will enable differen tiation between CDI 
(|Q'Neill. Beckwith fc Begelmanl [20121 ) and KHI contribu- 
tions (this work) to energy evolution in these systems. Fur- 
thermore, one could compare the non-linear evolution of tur- 
bulence examined here to similar turbulence that develops 
in joint KHI/CDI systems to determine how turbulent spec- 
tra, energy partitioning, and saturation levels differ between 
the two scenarios. 

The transfer function machinery developed and used 
here can be applied to other astrophysically relevant sys- 
tems. In particular, local simulations of magnetised accre- 
tion discs in the "mesoscale" regime (i.e., scales much larger 
than a ver tical scale height but much l e ss th an the disk 
radius) by ISimon. Beckwith fc Armitagei (|2012l) show that 
as larger disc scales are captured within the domain, tur- 
bulence driven by the magnetorotational instability (MRI; 
iBalbus fc Hawlevl [l998) develops structure on these larger 
scales at the expense of small scale structure. This behaviour 
is indicative of either an inverse cascade of energy or direct 
communication between small scales and large scales. In ei- 
ther case, applying our transfer function analysis to these 
mesoscale simulations will lead to a better understanding of 
energy flow in MRI turbulent disks. We are currently pur- 
suing this line of work. 

While astrophysical scenarios often lend themselves 
nicely to powerful computational studies, various obstacles 
(e.g., numerical convergence, multitude of important phys- 
ical processes, wide range in physical and temporal scales) 
force numericists to omit certain physics. When restricted 
to the ideal (i.e., inviscid and non-resistive) MHD limit, one 
often conjectures that numerical dissipation behaves suffi- 
ciently similarly to physical dissipation, even in situations 
where dissipation may be an important physical process for 
the problem at hand. For instance, dissipation of turbulence 
arising from the MRI is an important problem in compact 
object accretion disc physics, yet these studies are commonly 
performed in the ideal MHD limit. As an example, attempts 
to model the effect of the vertical dissipation profile on the 
emergent accretion disc spectrum are very i mportant for 
understanding observations of X-ray binaries |Turner|[200^ ; 
iHirose. Krolik fc Stonel l200rj ; iBlaes et all 120061 ). A reason- 
able question to ask is whether numerical dissipation leads to 
unwanted numerical artifacts in the absence of physical dis- 
sipation. Our work demonstrates that the details of numer- 
ical dissipation do not affect the physics of KHI-produced 
MHD turbulence on scales larger than the dissipation scale. 
Therefore, studies of ideal MHD turbulence conducted with 
codes comparable to Athena are not plagued by numerical 
effects due to the nature of numerical dissipation. 
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APPENDIX A: NUMERICAL ISSUES 

Here, we address some of the numerical issues in the KHI 
simulations with the Athena code in order to justify our 
particular setup. 



Al Riemann Solvers 

The available suite of Riemann solvers implemented in 
Athena for approximating physical fluxes across cell inter- 
faces are the Roe, HLLD, and HLLE solvers for magne- 
tohydrodynamics (MHD) and the exact, Roe, HLLC, and 
HLLE solvers fo r hydrodyn ami cs (for descrip tions of Rie- 
mann solvers, see lTorclll999l and lLevequekoOa ) . Waves trav- 
eling between grid cells are dispersive and wave propagation 
can be visualised as a Riemann 'fan', with the fastest-moving 
'leftward' and 'rightward' waves defining the fan edges and 
an ensemble of intermediate-speed waves composing the fan. 
The HLLD and HLLC solvers consider many intermediate- 
speed waves when computing fluxes, while the HLLE solver 
omits these waves and only accounts for the fastest waves 
propagating in each direction. Unlike the HLL- family of 
Riemann solvers, the Roe solver constructs the exact solu- 
tion to a linearised form of the equations at the cell inter- 
faces. 

Figure lAll shows the comparison of Riemann solver per- 
formance in Athena for the two-dimensional KHI in hydro- 
dynamics and MHD. While all solvers capture similar linear 
growth rates as measured from the magnetic energy evolu- 
tion, the HLLE solver diverges from Roe and HLLD dur- 
ing the onset of the non-linear evolution. Specifically, mag- 
netic energy in the HLLE run saturates at a level approx- 
imately 30% less than that of the Roe solver and spends 
a much longer time at a saturated state prior to entering 
the decay phase of the instability. Such diffusive behavior in 
the HLLE solver has been repor t ed in other contexts, too 



(Mignonc, Ugliano 


k Bodd 20091; Beckwith & Stone! 1201 ll: 


lO'Neill. Beckwith & 


; Begelman 2012)). all of which suggest 



that previous investigations of the magnetised KHI that rel; 
on the HLLE solver (e.g., iBucciantini 



ied rvm tnat rely 
Del Zannal boOrj ; 



IZhang. MacFadyen fe Wand 1200*91 ) may suffer from similar 
effects. In our KHI simulations, we use the HLLC and HLLD 
solvers exclusively for hydrodynamics and MHD, respec- 
tively, with the hllallwave configure option turned on to 
include the full interpolated Riemann fan. 
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Figure Al. Comparison of Riemann solver performance in 2D 
KHI simulations with resolution (N x X N y ) = (2048 X 2048) com- 
puted with the Athena code. Top panel: Time evolution of volume- 
averaged kinetic energy, (Er), relative to the volume-averaged 
total energy in the computational box, (-Etot)i for hydrodynamic 
simulations performed with the exact (green line), HLLC (black 
line), HLLE (red line), and Roe (blue line) Riemann solvers. The 
evolution of (Er) is essentially independent of the chosen solver 
in hydrodynamical simulations. Bottom panel: Time evolution of 
volume- averaged magnetic energy, (%), multiplied by a factor of 
5 (solid lines) and (Er) (dashed lines), each relative to (Btot), for 
MHD simulations performed with the HLLD (black line), HLLE 
(red line), and Roe (blue line) Riemann solvers. The saturation 
level of (Em) f° r the run adopting the HLLE solver is ~ 30% 
below that of the runs that used the HLLD and Roc solvers. 



A2 Linear Growth 

To demonstrate that our computational setup indeed pro- 
duces sensible linear growth of the KHI, we compared the 
development of a simple, equal density realisation of the 
KHI achiev ed with Athena to est imate s of linear growth 
provid ed in IChandrasekharl (|l96ll) and Miura fc Pritchettl 
(1982). The expression in IChandrasekharl 1)19611 ) describes 
the growth of an infinitesimally sharp shear layer in an in- 
compressible, weakly magnetised fluid as T ~ kUo (using 
our notation), which would correspond to a value of T ~ 3 
in our code units. The growth rates in iMiura fc Pritchettl 
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Figure A2. Volume-averaged rms velocity transverse to the 
shear layers for the resolved shearing model 3M512 (solid line), 
the discontinuous shearing model 3M512J (dotted line), and the 
extended domain models 3M512z2 (dashed line) and 3M512z4 
(dash- dot line). The linear growth rate for the discontinuous pro- 
file is dissimilar from those with resolved shear profiles. Simula- 
tions with an extended z-domain exhibit different behaviour from 
the smaller box fiducial run during the early stages of the non- 
linear decay phase of the instability. 



(1982) are more applicable to our setup in that they in- 
corporate a finite-width shear layer and compressibility, but 
also unfortunately rely on the approximation that the modes 
are short in wavelength compared to the box size, which is 
not satisfied for o ur k = 2tt/L perturbation s. The maximum 
growth rate from iMiura fc PritchettJ (|l982h most appropri- 
ate for our setup is V ~ 7, wh ich is considerably faster than 
that of IChandrasekharl (|l96ll ) because it occurs on a much 
smaller physical scale. Empirically, our fastest growth rates 
are measured to be T ~ 5, which falls comfortably between 
the two estimates. Furthermore, when we conducted KHI 
test cases featuring perturbations considerably smaller in 
scale th an L, we found growth ra tes more comparable to 
those in IMiura fc PritchettJ (|l982| ). We therefore conclude 
that the linear development of the KHI in our simulations is 
consistent with theoretical expectations for the instability. 



A3 Discontinuous vs. Resolved Shear Layers 

An inadequately resolved shear interface may result in the 
accumulation of numerical truncation error, causing unphys- 
ical realisations of the subsequent evolution. To quantify the 
degree to which the energetics are affected by the presence 
of an unresolved shear layer, we repeated the 3M512 simula- 
tion with the hyperbolic tangent interfaces for velocity and 
density replaced by jump discontinuities, 



Vy(z) 
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1, \ Z \ > Z 
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shear interfaces. Figure IA2] compares the time evolution of 
the volume-averaged root mean square (rms) velocity trans- 
verse to the shear layers, (u?) 1//2 , for the cases of resolved 
(solid line) and discontinuous (dotted line) interfaces. The 
initial onset of instability for 3M512J occurs sooner than 
in 3M512 because the accumulation of truncation errors at 
the interface permits perturbations at smaller scales (i.e., 
faster growth rates) than would be available for a finite- 
width shear layer. Despite the triggering of the KHI from an 
unresolved interface, the ultimate saturation and late-time 
evolution of 3M512J remains similar to that of 3M512. 



A4 Extending the Domain Transverse to the 
Shear Layer 

The choice of periodic boundary conditions was motivated 
by its ease of implementation and its attractive consequence 
of energy conservation within the domain. As the KHI 
evolves to the non-linear regime, propagating waves and 
fluid that exit through one boundary will re-enter through 
the opposite boundary and interact with the flow. A poten- 
tial concern is that cross-boundary interactions may sub- 
stantially affect the evolution of the flow and produce an out- 
come driven by numerical, rather than physical, processes. 

Figure IA2I compares the evolution of the volume- 
averaged rms velocity transverse to the shear layers in model 
3M512 with those of models 3M512z2 and 3M512z4, for 
which the z-domain is extended by a factor of two and four, 
respectively. For the 3M512z4 simulation, the z-domain is so 
far extended that the turbulent regions are isolated in z and 
the turbulence never crosses the z-boundaries. The data of 
Figure [A2l show that extending the domain transverse to the 
shear layer does not substantially affect the linear growth or 
the peak (v 2 ) 1 ^ 2 amplitude that is achieved. A difference 
arises only when the KHI is well into the non-linear decay 
phase. At this point, the flows are both very turbulent, but 
(v 2 ) 1 / 2 in the smallest domain decreases more rapidly after 
the peak amplitude than it does in the larger domains. Af- 
ter this, however, the two decay rates become approximately 
parallel, suggesting that cross-boundary interactions do not 
grossly affect the asymptotic shape of the decay phase even 
if they do adjust its levels. 



APPENDIX B: DERIVATION OF SPECTRAL 
ENERGY TRANSFER FUNCTIONS 

Spectral energy transfer analysis was fi rst intro- 
duced in the incompressible limit by iKraichnanl 
(1967). Transfer analysis is a well-developed tool 
for studying MHD turbulence in the incompressible 
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2005) and compressible (|Fromang & Papaloizou 




2003 


Fromang et al. 20071; Simon. Hawlev & Beckwith 




2009|; 


Pietarila Graham, Cameron & Schiisslerl 2010f) 




lim- 



We refer to this simulation as 3M512J, where the J refers to 
the jump discontinuities in velocity and density across the 



its. Transfer th eory was outlined for compressible 
MHD in IPietarila Graham. Cameron fc Schiisslerl 
(2010), which was a generalisation of the incom- 

J ressi ble treatment of lAlexakis. Mininni fc Pouquej 
2005) . Here, we expand on the transfer analysis of 
IPietarila Graham. Cameron fc Schiisslerl (|2010l ) by incorpo- 
rating a decomposed velocity; thus, separating the transfer 
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mechanisms involving the velocity field into components 
due to the background shear flow and turbulent motions. 
This allows one to distinguish between energy transfer 
arising due to turbulence from that due to the background 
flow. 

The basic philosophy behind deriving the transfer func- 
tions is to start by taking the complex conjugate of the 
Fourier transform of the conservation equations to ob- 
tain time-evolution equations of energy densities in Fourier 
space. The Fourier transformed conservation equations are 
then dotted with the Fourier transform of the appropriate 
quantity. The result is the time derivative of a spectral en- 
ergy density being equated to many individual terms. These 
terms are the transfer functions and describe energy transfer 
from one energy reservoir to another, mediated by a force. In 
what follows, we derive the magnetic, kinetic, and internal 
energy transfer functions, each in turn. 

The primitive form of the induction equation is, 



dB „ . 
-=Vx(vxB), 



(Bl) 



where B is the magnetic field and v is the fluid velocity field. 
We decompose the velocity field into a turbulent velocity, vt, 
and a shear velocity, v sh , according to, 



V = V sh + Vt, 



where, 



v sh = «sh(z)y 



v y (x,y,z)dxdy. 



(B2) 



(B3) 



Replacing the velocity field in Eg uation IB 1 1 with the decom- 
posed velocity defined by Equation IB2I taking the complex 
conjugate of the Fourier transform, where the Fourier trans- 
forrrlj of a quantity /(x) is given by, 



F(k) 



/(x)e 



-zk-x i3 

a x, 



(B4) 



and dotting the result with B(fc), we obtain the equation 
representing the transfer of magnetic energy in fc-space, 

dE M (k) 



dt 



Tbba(A;) + 5*BBA(fc) + Tkbt(A;) + Skbt(A;)- 



T KB p(fc) + D M (k). 



(B5) 



The left hand side is the time-derivative of the spectral mag- 
netic energy density, where, 



E M (k) = -B(fc)-B*(fc). 



(B6) 



The terms on the right-hand side of Equation IB5I are 
the magnetic energy transfer functions. We first de- 
scribe the transfer function notation and then identify 
each term explicitly. Transfer functions with the nota- 
tion Txyf(&) depend only on turbulent velocities, v t , 
those with notation Sxyf(&) depend only on the back- 
ground shear velocity, v s h, and those with notation 
Xx_YF(k) have a mixed velocity dependen c e. As described by 
IPietarila Graham. Cameron fc Schiisslerl (|20T0t ). the trans- 
fer function T, S, Xxyf(&) measures the net energy trans- 
fer rate from all scales of reservoir X to scale k of reservoir 



6 Fourier transforms of a quantity are denoted by a 'hat', , 
not to be confused with those hats implying unit vectors (e.g., 
x). 



Y, where the energy exchange is mediated by the force F. 
The net energy transfer from reservoir X into reservoir Y at 
scale k is positive (negative) for T, S, XxYF(fc) > (< 0). 
In other words, energy is lost by reservoir X and gained by 
reservoir Y at scale k for T, S, Xxyf(&) > and vice versa 
for T, S, XxYF{k) < 0. The available energy reservoirs are 
kinetic (K), magnetic (M), and internal (I). The cascade of 
magnetic energy to other scales from within the magnetic 
energy reservoir is described by the terms, 

T BBA (k) =-B(fc)-[^Tv)Bi*(Jfe) (B7) 

S B BA(fe) = -B(fc)-[(vI h ^V)B]*(fc). (B8) 

The transfer of energy from the kinetic energy reservoir to 
the magnetic energy reservoir by turbulent (T) and shearing 
(S) fluid motions twisting and stretching field lines are, 



T KBT (fc)= B(fc) ■ [(B • V) v t ] (k) 



(B9) 



S'KBT(fc) = B(fc) ■ [(B ■ V) v sh ] (k). (BIO) 

Magnetic energy transfer from the kinetic energy reservoir 
by compressive motions via the magnetic pressure force is 
given by, 

T KB p(fc) =-§(*) -jBCV^OH*), (Bll) 

where there is no 5*KBp(fc) transfer function with a back- 
ground shear velocity dependence because, 

dv s h{z) 



V sh = 



dy 



= 0. 



(B12) 



Formally, the magnetic energy transfer function expression 
(Equation lB5[) is analytically exact in the omission of the nu- 
merical magnetic energy dissipation term, DM(k). However, 
numerical schemes have dissipative effects. Therefore, any 
inequality that arises from the sum of the transfer function 
terms on the right-hand side is folded into Divi(fc), which is 
a measure of the numerical magnetic dissipation. The other 
transfer function equations will have associated numerical 
dissipation terms as well. 

The kinetic energy transfer functions are derived in a 
similar fashion as was done for the magnetic energy transfer 
functions. Starting from the primitive form of the momen- 
tum equation, 
<9v 



-p (v • V) v — VP + (V x B) x B, (B13) 
and the conservative form of the momentum equation, 

d( P v) 



dt 



= -V 



pvv 



BB 



(B14) 



the velocities in these equations are decomposed according 
to Equation lB2l Here, p represents the mass density, P is the 
pressure, and I is the identity matrix. The complex conju- 
gate of the Fourier transform of Equation lB13l is dotted with 
pv and the complex conjugate of the Fourier transform of 
Eouation lB14l is dotted with v. Combining these two result- 
ing equations yields the expression representing the transfer 
of kinetic energy in fc-space, 

dE K (k) _ TlKC ( k ) + SlKC ( k ) + TKKA ( k } + x KKA (k) + 

TBKT(fc) + SBKT(k) + T B Kp{k) + SBKp(fe) + 

TkkcW + S K Kc(k) + X KK c{k) + D K {k), 

(B15) 



dt 
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where the spectral kinetic energy density is defined by, 

E K (k) = ~ (v(fc) ■ [H*W + IH(fc) ■ v*(*)) • (B16) 

The transfer functions describing the exchange of kinetic en- 
ergy from within the kinetic energy reservoir by compressible 
motions due to turbulence and background shear are, 



Tkkc (fc) 



\ (*(*) 



[v t (V-pv t )] (fc) 



(B17) 



S K Kc(fc) = - i (v£(fc) ■ [vT(V~^)]*(fc)) , (B18) 

respectively. The corresponding cross-term transfer function 
is, 

^KKCO = 

- \ (jt(k) ■ [v sh (V ■ pv t )] (fc) + v^(fc) • [v t (V ■ pv t )] (fc) 

1 
2 
1 
2 



The transfer of energy within the kinetic energy reser- 
voir by advection is described by the transfer functions, 

TkkaW = 
i 



v sh (*0 • Kh (V ■ pv t )] (fc) + vt(fc) ■ [v t (V ■ pv sh )] (fc) J 
vt(k) ■ [v sh (V ■ pv sh )] (fc) + C^(fc) • [v t (V ■ pv sh )] (fc) ) . (B19) 



[pvt] ■ [(v t ■ V) v t ] (fc) + vt(fc) ■ [(pvt ■ V) v t ] (fc) 



(B20) 



x KKA( k ) = 
1 
2 
1 



2 

i( 

1 
2 



■ [(v t 


"v)~v7]* 


[(Vsh 


■ V) v t ] 


■ [(v sh ■ V) v t ] 


• [(v t 


■ V) v Bh ] 


[(vt ■ 


V) v Bh ] 



\ UPV^JW ' [(Vsh • V)v t ] (fc) + Vsh(fc) • [(v s h-V)v t ] (fc)) 



[pv t ](fc) ■ [(v t • V) v Bh ] (fc) +v t (fc) ■ [(pv t ■ V)v ah ] (fc) . (B21) 



Generally speaking, the transfer functions involving 
mixed velocity terms, which are denoted by XxYF(fc), are 
not intuitively graspable. However, these cross-terms only 



appear for the transfer functions describing the energy cas- 
cade within the kinetic energy reservoir by compressive mo- 
tions and advection. Energy transferred from the magnetic 



energy reservoir to the kinetic energy reservoir by fluid mo- 
tions via the magnetic tension force are, 



TbktCO 



i(bvt]( fc) .[I 

~ ■ [~ (B ■ V 



(B-V)B (fc) + vt(fc) • [(B ■ V) B] (fc) 



(B22) 



)B| (fc) + v sh ■ [(B ■ V)B] (fc) 

(B23) 

and by compressive turbulent fluid motions via the mag- 
netic pressure force are, 

TbkpCO = - | f [pvj](fc) ■ [^VflsJ ( k )+w t ■ [^VB 2 j (fc)J (B24) 

s BK p(fc) = - i (jp^iW ■ [^ vs2 ] + ^ ■ [~ vs2 

(B25) 

Finally, the transfer functions describing energy ex- 
change from the internal energy reservoir to the kinetic en- 



ergy reservoir by compressive motions are, 

IteoW = - i (jpVt\(k) ■ Jj VP ] (fc) + v^ ■ [VF]*(fc)^ (B26) 
SiKc(fc) = - i ^[pv^j(fc) • [^ VP ] (*)+VA-[VP]*(fe)J. (B27) 



Again, following the same procedure as for deriving the 
magnetic and kinetic energy transfer functions, we start with 
the internal energy equation, 



^_ v . V P + 7 P(V-v), 



(B28) 



where 7 is the adiabatic index. Decomposing the velocity 
in this equation according to Equation IB2I and multiplying 
the complex conjugate of the Fourier transform by P(k), the 
equation describing the transfer of internal energy in fc-space 
becomes, 



dEi(k) 
It 



= TkiaO) + SkiA(fc) + TkicO) + Di(k). (B29) 



The internal energy density in fc-space is defined as, 

P(k) 



7' 



(B30) 



and the transfer functions associated with energy exchange 
from the kinetic energy reservoir to the internal energy reser- 
voir by advection and compressive motions are, 



T K1A (k) =• 
SklA(fc) =• 
T K ic(fc) =- 



1 



7-1 
1 



7" 



VP(k) 
VP(k) 



1 

VP 



(v t ■ V) P 



(v sh • V) P 



(fc) (B31) 



7' 



■VP(k) 



-^P(V.vO 



(fc). 



(B32) 
(B33) 



The transfer function analysis presented above can be 
extended to the case of a viscous and resistive fluid. We 
derive these additional dissipation transfer function terms 
follow ing t he procedure outlined in iFromang fc Papaloizoul 
(|2007l ) and lSimon. Hawlev fc Beckwithl (j2009h . Turning our 
attention first to the induction equation, the inclusion of 
Ohmic resistivity introduces the term, 7?V 2 B, to the right- 
hand side of Equation IB1I where 77 is the resistivity. Tak- 
ing the complex conjugate of the Fourier transform for this 
Ohmic dissipation term and dotting it with B(fc) yields, 

T„(fc) = V (§(fc) ■ [V2B]*(fc)) (B34) 

Incorporating viscosity would add the terms, (V • r) / p, 
and, (V • t), to the right-hand sides of Equations IB13I and 
IB14I respectively, where the stress tensor for an isotropic 
fluid is, 



nj = 2/i (eij - - (V- v)^-) , 



(B35) 



where fi is the dynamic viscosity, Sij is the Kronecker delta 
function, and the strain rate tensor is given by, 



e« = \ [(Vv) 



(Vv) 



(B36) 
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One can show that, 




(B37) 



The dynamic viscosity is related to the kinematic viscosity, 
v, by v = fi/p. Taking the complex conjugate of the Fourier 
transform of the viscous term and dotting the conservative 
form with v(fc) and the primitive form with [pv](fc) gives 
the transfer function describing the viscous dissipation, 



Note that due to the non-linear nature of the resistive 
and viscous terms, we choose not to decompose the velocity 
field in the definitions of T v (k) and T v {k) in order to simplify 
their interpretations. 




(B38) 



© 0000 RAS, MNRAS 000, 000-000 



